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Abstract 


We have derived stability results for high-order finite difference approximations of mixed 
hyperbolic-parabolic initial-boundary value problems (IBVP). The results are obtained 
using summation by parts and a new way of representing general linear boundary con- 
ditions as an orthogonal projection. By rearranging the analytic equations slightly we 
can prove strict stability for hyperbolic-parabolic IBVP. Furthermore, we generalize our 
technique so as to yield strict stability on curvilinear non-smooth domains in two space 
dimensions. Finally, we show how to incorporate inhomogeneous boundary data while re- 
taining strict stability. Using the same procedure one can prove strict stability in higher 
dimensions as well. 




1 Introduction 


When solving a partial differential equation numerically it is necessary to have some bound 
of the growth rate of the solution, since otherwise round-off errors could grow arbitrarily 
fast. This upper bound can be established by ensuring some kind of stability. We have 
elected to use the energy method, because it can be applied to the continuous as well as 
the discrete model. Furthermore, it can be applied to general domains, which is important 
when studying multidimensional problems. 

Stability of the continuous problem is established by means of an integration-by-parts 
procedure introducing boundary terms, some of which must be eliminated to ensure sta- 
bility. For the finite difference model integration by parts is replaced by summation by 
parts. This amounts to designing the discrete difference operator ensuring that, in ad- 
dition to the accuracy requirements, certain conditions of antisymmetry are met. As a 
consequence, the common problem of finding proper ” numerical” boundary conditions 
will be eliminated; they will be built in the discrete difference operator. 

The analytic boundary conditions are yet to be incorporated. We propose a certain 
projection operator, which interacts with the difference operator so as to generate bound- 
ary terms that are completely analogous to those of the continuous problem. This can be 
done for any type of linear boundary conditions. Thus, an energy estimate is obtained 
for the discrete problem, provided there is one for the analytic model. This conclusion 
remains true for domains in several space dimensions, even if the boundary is non-smooth. 
Furthermore, using this projection operator allows us to derive stability results for a larger 
class of finite difference operators than those considered in [5]. Stability will be proved 
for high-order finite difference approximations of mixed hyperbolic-parabolic variable co- 
efficient systems subject to general inhomogeneous boundary conditions. 


1.1 An Introductory Example 

To illustrate the underlying principles of the energy method we consider the convection- 
diffusion equation 

lit — u X x + u x , x 6 (0, 1 ) t > 0 
u(a:,0) = f(x) 
u(0, t) = 0 

«*(M) = 9(t) 

In the sequel we shall use the standard L 2 -scalar product 

( u , v)= uvdx 

Jo 

with the corresponding norm defined as ||u|| 2 = (u,u). 

We can obtain an a priori estimate for this example using the following tools. 
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(i) Integration by parts: 


— ||u|| 2 = 2 (u,u xx ) + 2(u, u x ) = — 2||u a: || 2 + 2(ti, u x ) + 2uu x \l 
at 


(ii) Boundary conditions: 


-||u|| 2 = -2|K|| 2 + 2(u,u,) + 2u(l,t)g(l,t) 
at 


(iii) Cauchy-Schwarz inequality: 


d 


J t \ lull 2 <'-2||u I || J + 2||u[||lu x || 4 2 u(l,()<l(M) 


(iv) Algebraic inequality: 
implies (e = 1) 


2\xy\ < ex 2 + e 1 y 2 


j t \\u\\ 2 <-\\u x \\ 2 + \\u\\ 2 + u(ht) 2 + g(ht) 2 

(v) Sobolev inequality: 

K<<IM 2 + (e-' + l)IM| 2 

is used to eliminate u(l,t) (e = 1) 

d 


dt 


IMI 2 <3||«|| 2 + <Km) 2 


which can be solved analytically to yield 

IW-,i)H 2 <c 3 ‘ (ll/ll 2 + J\(r) 2 dr) 

If we are to obtain such an estimate for a system of equations we will also need 


(vi) The adjoint of A: 


(u, Av) = (A r u, v) 


Summing up, the energy method boils down to the six basic “tools” above. In the subse- 
quent sections we shall see how these principles can be modified so as to give an energy 
estimate for the semi-discrete system. 
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2 General Principles for the Semi-discrete Case 


In this section the basic principles of the energy method will be transferred to the semi- 
discrete case. Furthermore, a number of lemmas, which will be needed later, will be 
proved. Throughout this section grid vectors will be denoted by v T = (vq . . . uj), v 3 6 R rf . 
Difference operators approximating d/dx will be designated by 


/ d 00 I 



\ d u $I 


dovl ^ 

dw I ) 


I G R 


dxd 


where D is written as a square matrix for convenience; in reality D will be a banded 
matrix, where the bandwidth is independent of the mesh size h = Ljv, 


2.1 Summation by Parts 

In the semi-discrete case we employ summation by parts instead of integration by parts. 
The basic idea is to use difference operators satisfying 

(u,Dv)h = ulvv - ulv 0 - (Du,v)h (1) 

with respect to a weighted scalar product 


{u,v)h - hY^ VijU-Jvj 

i,j=0 


It should be remarked that the usual Euclidean scalar product cannot be used. To prove 
the existence of summation by parts, it suffices to consider scalar products on the form 


S = 


E' 1 ) \ 

I 

/ 


E (,) € R( r,+1 ) dx ( r ' +1 )' f j / = 1,2 


( 2 ) 


where the blocks of E are given by = cr^/, / € R dxd ; n and the elements of E^, 
l = 1,2 are independent of h. The following existence proof can be found in [5]. 


Proposition 2.1 There exist scalar products (2) and difference operators D of accuracy 
2p — 1 at the boundaries and 2 p in the interior, p > 0, such that the summation-by-parts 
property (1) holds. 


Confining ourselves to the case where E^ and E^ 2 ^ are diagonal we have the following 
existence theorem [4]. 
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Proposition 2.2 There exist diagonal scalar products (2) and difference operators D of 
accuracy p at the boundaries and 2 p in the interior, 1 < p < 4, such that the summation- 
by-parts property (1) holds. 


Remark: If one omits the requirement that the boundary stencils be at least accurate 
of order p for a given interior accuracy 2 p, it is possible to prove summation by parts for 
diagonal scalar products and difference operators D of arbitrary order of accuracy [7]. For 
a given boundary accuracy p , however, it may be necessary to resort to interior stencils 
of accuracy q 2 p, which may render these operators useless in practice. 

The actual computation of the operators above is ill-conditioned, since it involves the 
solution of a rank-deficient problem. Using a symbolic language it is possible to solve 
for D exactly, the elements of which in general will depend on one or more parameters. 
Explicit examples can be found in [6]. For details on the algorithms we refer to [8]. The 
simplest example is furnished by 

/ -i i 
-0.5 0 0.5 




V 


-0.5 0 0.5 

-1 1 


(3) 


with the corresponding scalar product 

( 0.5 

£ = h 


0.5 


(4) 


Summation by parts can be generalized to several space dimensions if we restrict our- 
selves to diagonal norms. To simplify the notation we consider only the two-dimensional 
case. A general proof is given in [6]. The grid function is partitioned as u T — 

( U T...ul), U J = ( 


U t 


0 j 


, u 


VM 


), j = 0, . . . , v 2 . Define the weighted scalar product as 


(u,v) h = hJ2Y, (T i (T 3 U Ij V ' 

t=0 j= 0 




(5) 


where we have assumed the same number of grid points in both dimensions for conve- 
nience only; h = hih 2 is the cell area. Let Di and D 2 denote the difference operators 
approximating d/dx i and d/dx 2 . Define 


1 * \ » 

{D\u)ij = t~Y1 ( D 2 u)ij = — ^2 djkUik 


Ais ; 


(6) 


‘2 k = 0 
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where it is assumed that the cr’s and d' s satisfy (1). Hence 

V / V V 

(«, Div) h = h 2 Y a j ( X (TlU ^j X d ik y k] 

j - 0 \i=Q k= 0 

and a similar expression holds for (u, The parenthetical expression satisfies (1) for 

each j. We thus arrive at 


Proposition 2.3 Let the discrete difference operators D\ and D 2 be defined by (6). Sum- 
mation by parts then holds in both dimensions 

V V 

(u, D\v) h = h 2 X ~ h 2 X a J u oj v oj - (Am, v) h 

j — 0 j = 0 

1 / u 

(u, D 2 v) h = h x Y VtuJvViv - hi Y a i u Io v io - ( D 2 U, v) h 

1=0 i = 0 


where (v)h defined by (5). 

Remark: This is the discrete counterpart of the two-dimensional divergence theorem. 
With a general domain f 1 we assume that there is a smooth map £ = £(x) taking LI 
onto the unit cube where proposition 2.3 can be applied. The assumption of such a 
map £ is necessary in order for finite difference methods to apply to curvilinear domains. 
Consequently, integration by parts can always be replaced with summation by parts in 
the discrete case. It is presently unknown if it possible to obtain the summation-by-parts 
property in more than one dimension using non-diagonal norms. 


2.2 Projections 

Suppose that the model equation of section 1.1 were discretized as 

v t = D 2 v + Dv , 

v(0) = f 

where we have assumed homogeneous Neumann data for convenience; it will be shown 
later how to treat inhomogeneous boundary conditions. For every fixed h the problem 
above is a constant coefficient ODE system with a unique analytic solution. Consequently, 
there is little hope that the discretized boundary conditions vo(t) = (Z)u) l/ (t) = 0 are 
fulfilled, since they have not been accounted for so far. 

Denote by V C Ft"" 1 " 1 the vector space where Uo(0 = {Dv) v (t) = 0, and let P be a 
projection of v onto V. Multiplying (7) by P yields 

(Pv) t = P ( D 2 v + Dv) 


5 



Any solution satisfying the boundary conditions must obey v = Pv, whence 

v t = P {D 2 v + Dv'j 

Conversely, we have 

Proposition 2.4 Let P £ 7? sXs be a given projection independent of t, and suppose that 
v(t) £ R s is a solution of the non-linear ODE system 

v t = PR(t,v) + (7 - P)g t 
v(0 ) = f 

where f satisfies f = Pf + (7 — P)g(0). Then 

v(t) = Pv(t) + (7 - P)g(t), t > 0 


Proof: 

Since P is independent of t , premultiplication of (8) gives (P 2 = P) 

(Pv) t = PR(t,v) 


Using this equality in (8) implies 


v t = (Pv + (I - P)g) t 


Hence, by integration, 

(I - P)(v(t) - g(t)) = (I - P)(f - g(0)) 

which proves the proposition. O 

Remark: g(t ) represents the boundary data, and (7 — P)(v — g) — 0 is the extension of 
(7 — P)v = 0 to inhomogeneous boundary data. Proposition 2.4 thus tells us that any 
solution to (8) will satisfy the boundary conditions if the initial data do so. 

In general P is not uniquely defined. Consider the vector space V = {u £ R I/+1 |i>o = 
0,1V = Vv-i}. Then 



0 

1 

\ 


0 

1 

\ 

p = 


1 

1 0 ) 

p = 


0 1 
1 / 


both imply Pv £ V. To shed some light on how to choose P, we apply the energy method 
to (7) 

j t \\v\\l = 2(v,P(D 2 v + Dv)) h 
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If P were self-adjoint w r t (•,•)&, then 

^||v||fc = 2(Pv,D 2 v + Dv) k = 2(v,D 2 v + Dv) h 

where the last equality follows from proposition 2.4. The crucial condition to obtain this 
equality is expressed by 

(u,Pv) h = (Pu,v) h (9) 

which states that P is an orthogonal projection (using the weighted scalar product (•, -)h). 

Suppose that u(x,t ) € R rf , x G R" is a solution to 

u t — F(x, t, d)u x € fl 
L(x , d)u = 0 x €E r 

where d denotes the n-dimensional gradient; T is the boundary of fl. This system is 
discretized in space, possibly requiring a coordinate mapping onto the unit cube 

v t = PG(t, D)v 

The projection P should be such that v fulfills 

L t v = 0 

where L now represents a discretization of the analytic boundary conditions. Let V — 
{v € R m |L r t; = 0}. According to the preceding discussion P is taken to be the orthogonal 
projection onto V (with respect to ■)*.)- The boundary conditions can be written as 

Q t T,v = 0 

where Q = Hence, the boundary conditions are fulfilled for all vectors v that 

are orthogonal to the column space of Q, the orthogonal projection onto which reads 
Q(Q t TiQ)~ 1 Q t 'L. In case E = I this is the standard projection. The desired boundary 
projection is thus given by 

p = /-Q(g T Eg)- 1 g T s 

or 

P = /- E“ 1 Z,(L t E'U)- 1 L t (10) 

Remark: In order for the projection to be well-defined the inverse of L T E _1 L must exist, 
which follows iff L has full column rank. The latter will follow from assumptions on the 
analytic boundary conditions (consistency arguments). 

Proposition 2.5 Suppose that L has full column rank, and let P be defined by (10). 
Then 

(i) P 2 = P 
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(ii) £ P = P T E 

(iii) v - Pv <£=$■ L? v = 0 


Proof: 

All statements are immediate consequences of (10). 

Remark: The second statement of proposition 2.5 is equivalent to (9). 


2.3 A Discrete Sobolev Inequality 

As seen in section 1.1 it is necessary to have a Sobolev inequality. The following proposi- 
tion shows that there is a discrete Sobolev inequality for the norms that we are interested 
in. We present it in a form suitable for proving strict stability. 

Proposition 2.6 Let || • ||^ and D be defined by (2) and (1), respectively. Then 

ML < ell-Di’lli; + («'* + 1 + O(h)) IMIJ 


where e > 0. 

Proof: 

Choose k , l such that 

\v k \ 2 = min,(h| 2 ) 

M 2 = max, (h| 2 ) = |u|^ 

Eq. (2) implies that 

IMIft > fc(Ai|u {1) | 2 + A 2 |u (2) | 2 ) +/i J2 hi 2 

j=n+l 

where A 1)2 >0 are the smallest eigenvalues of X^ 1 ’ 2 ). Note that A 1>2 are independent of h. 
Hence 

IM|^>(l-ft(r 1 (l-A a ) + r 2 (l-A 2 )))h| 2 

where we have used hi/ = L = 1. If c = r x (l - A x ) + r 2 (l - A 2 ) < 0 one immediately gets 
hi 2 < ||t>||j;. Otherwise we choose h such that he < 1. Hence 

+ K = rri ^ (11) 

for h < ho, where h 0 is a fixed number such that h 0 c < 1. 
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Next, we define a family of norms, which is obtained by shrinking the interior of (2); 
remain constant. Allowing a slight abuse of notation we write these norms as 

S 

(U,v)h,r,s = EWTVJ 

J=r 

where r > 0 and s < v. Shrinking the interior of D accordingly one has 

(u, Dv)h,k,i = N 2 - |u/t | 2 - ( Dv , v)h,k,i 


1 . e., 


ML ^ K| 2 + 2||£)u||fc,fc,i||v|U,fc,/ 


Obviously ||u|| fcifci / < || v || a , 0 ,„ = |MU, whence 

where (11) and the standard algebraic inequality have been used. 


□ 


2.4 Adjoint Operators 


As usual A T denotes the transpose of A. We know from section 1.1 that (u,Av) = 
( A T u,v ), i. e., the transpose of A is the adjoint operator. The question is whether A T 
also is the adjoint operator with respect to as defined by (2). Let 


A = 


^ A 0 



Aj = A(jh), j = hv = 1 


( 12 ) 


denote the matrix representation of A(x ) G x G [0, 1], Smoothness will be assumed 

as needed. 


Lemma 2.1 Let E and A be defined by (2) and (12), respectively. Then 


I(u,y4u)fc - (A T U,v)h\ < <3(A)|MU|MU 


Proof: 

Denote the commutator of E and A by [E, A], Then 

(it, Av)h = ( A t u , v)h + /m T [E, A]v 


where 


[S, A] = 


( pi'U"'] 


[EW./ll 2 '] 


( 13 ) 
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with 


[£(!), y^)] = 


001 (Ax — Ao) • • • a O r(A r — Ao) \ 
0 ... ai r (A r — Ax) I 


/ o 

—CTq\(A\ — Aq) 


\ — cr 0r (i4 r — Ao) — (J\r{A r — /4i) 


The other non-zero block has a similar structure. Assuming that A(:r) is differentiable we 
can apply the mean value theorem 


/ 


[E (1) ,A (1) ] = h 


0 

-(t 01 A[ 0 


ctoiA'k) ... cr 0r rA' r0 \ 

0 ... a lT {r-\)A' rA 


\ -<jQ T rA' r o -<Tir(r ~ 1 )A' rX 


0 I 


Hence, 

t(u, Av ) h - (A T u,v) h \ < c|A'UA 2 (|k <1) ]|o ( 1) I + |u (!, l|f ,!, l) < 0(*)|MUI|®IU 


which proves the lemma. □ 

Remark: According to lemma 2.1 the transpose of A is an approximate adjoint with 
respect to (•,•)/,, and the perturbation consists of lower order terms. The following as- 
sumption will be crucial when proving strict stability for hyperbolic systems. 


Assumption 2.1 Let A and S be given by (12) and (2). Then one of the conditions 
below is assumed to hold. 

(i) E is diagonal diag(cr 0 / . . . cr u I). 

(ii) The blocks of A satisfy Ao = ... = A rj and A l/ - T2 +i = . . • = A„. 

Corollary 2.1 If assumption 2.1 holds , then (u,Av)h = ( A T u,v)h . 

Proof: 

In both cases we get [E, A] = 0. The result follows immediately from (13). □ 

Remark: The latter criterion is satisfied if there is a 6 > 0 such that A(x) = const for 
0 < x < 6 and 1 — 6 < x < 1, and if h is chosen such that hr < 6, where r = max(ri, r 2 ). 


Corollary 2.2 Let A be given by (12). If A is symmetric then then ( u,Au)h = (Au,u)h- 

Proof: 

A symmetric implies that [E, A] is antisymmetric. Hence u r [E, A]u = 0. □ 
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2.5 Some Operator Estimates 


In this subsection we have gathered some operator estimates that will be needed in sub- 
sequent sections. The results are valid for norms defined by (2) unless otherwise stated. 
In particular, the estimates will be given in a form suitable for proving strict stability of 
the semi-discrete systems. 


Lemma 2.2 Let £ and A be defined by (2) and (12). Then 

\(u,Av) h \ < |>4|oo (1 + 0(h)) IHlfclMIfc 


where |v4|oo = sup |j4(:e)|. 


Proof: 

The definition of (-,•)* implies that (u,Av)h — hu T Av , where u = £ ^ 2 u, v = S ] / 2 u, and 
A = E 1 / 2 AE _1/ ' 2 . Taylor expansion yields A — A + R, 


( RW 


R 


RW ) 


with RW = 0{h), 1 = 1,2. Thus 

l(«, A«)»| < |/t|„||u||||ii|| + 0(h) (||i' 1 l||||6' 1 >|| + ||i' 2 >||||5' 2 »||) 


i. e., 


|(u,i4t;)fc| < (|A|oo + O(h)) ||u||||u|| 


where || • || denotes the standard Euclidean norm. Since ll«ll = IM \h, ||v 
lemma follows. 


\\v\\h, the 
□ 


Corollary 2.3 If, in addition to the hypotheses of lemma 2.2 , assumption 2.1 is fulfilled , 
then 


|(u, Av) h \ < lAUlHUHHU 


Proof: 

The hypotheses imply that A = A, and the corollary follows. □ 

Remark: Lemma 2.2 states that the growth rate induced by low order terms is the same 
(modulus L?(A)-terms) in the continuous and the semi-discrete case. 

It is well-known that (u, [D, A]u)^ < ||[Z), J 4]||/ l ||u||/ l ||t;||/ l , where ||[L>, A]||/i can be 
bounded independent of h. This result can be sharpened under certain circumstances. 
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Lemma 2.3 Let D be a difference approximation satisfying the summation-by-parts rule 
(1) with respect to a weighted norm (2), and define A by (12). Suppose that assumption 
2.1 holds. If A is symmetric, then 

(u,[D,A]v) h < p{[D, j4])||tt||/»||f IU 


Proof: 

According to the definition of the operator norm we have 

\\[D,A]\\l = max \\[D,A]v\\l = max hw T C T Cw 
IMIh=i IMI=i 

where w = E J / 2 u, C = Y}^ 2 [D, A] E -1 / 2 : Because of the assumptions on A (or E) we have 
C = DA — AD, where D — E^-DE -1 / 2 . Summation by parts implies that 


-/ 


E D = D s + D a , D s — — 


/ € R 


dxd 


and D a is an anti-symmetric matrix. Consequently, 

C = [E- 1/2 £> a E~ 1/2 ,A] 

where we have used [E -1 / 2 D S E -1 / 2 , A] = 0. Since D a is anti-symmetric and A symmetric 
we have C T = C , i. e., 

\\[ D , A ]\\h = mfx hw T C 2 w = p(C ) 2 

Finally, C = E 1 / 2 [D, A]E -1 / 2 implies that 

\\[D,A)\\ h = P ([D,A}) 


which proves the lemma. 


3 Homogeneous Boundary Conditions in One Di- 
mension 

We shall successively consider hyperbolic, parabolic and mixed hyperbolic-parabolic sys- 
tems. Variable coefficient matrices will be allowed. To simplify the presentation we shall 
only deal with the lower boundary x = 0, which is justified if we take the solution to have 
compact support. In general, the upper boundary x = 1 is treated in a fashion similar to 
the procedure at the lower boundary. 
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3.1 Hyperbolic Systems 

Consider the hyperbolic system 

u t = A u x + Bu + F x £ (0, 1) 
u(z,0) = f(x) A (x,t) = 

u_(0,<) = Lu + (0,<) LGR rflXd2 

where u £ R^, d] + <f 2 = d; A_, A+ is the partitioning of A into negative and positive 
eigenvalues. It is assumed that the elements of the diagonal matrix A never change sign 
at the boundaries x = 0 and x = 1 , and that there is a constant 7 > 0 such that 
A f) < —7 and A+(j,<) > 7 , j = 0, 1 This implies that the rank of L is constant. 
Furthermore, L is assumed to be “small”. 

The discrete boundary conditions are written as L T v — 0, where 

L t =[lI 0 ... 0 ) E R dlX(l/+1 ) rf (15) 

Here L£ = ( I —L ) € Rr 1 , the latter L being the analytic boundary operator. It 
follows immediately that rank(T) = rank(7) = d\. The hypothesis of proposition 2.5 is 
thus satisfied, and we have the semi-discrete system 

v t = P(ADv + Bv + F) 
v{0) = / 


A(0,f) 


A = 


A (1,0 ) 


(16) 


A 


A +(x,t) 


( 14 ) 


Proposition 3.1 Let (•,•)/, be given by (2) and suppose that D satisfies the conclusion 
of proposition 2.1. If P is defined by (10) and (15), then the solution of (16) satisfies an 
energy estimate 

HOIK + £ (Mr)l 2 + Mr)| 2 ) dr < A'e<"' +£ W>' (||/||J + £ ||F(r)|| ! t dr) 

Proof: 

The energy method yields (using propositions 2.5, 2.4) 

4lMli[ = 2 i v , v t)h = 2(w, P(ADv + Bv + F)) h = 2(n, A Dv) h + 2(v, Bv) h + 2 {v,F) h 
at 

Summation by parts implies ( v v = 0) 

(u, A Dv) h = -ujAono - ( Dv , Av) h - (u, [ D , A]v) h 
Hence, by lemma 2.1 

(v,ADv) h < -^a 0 u 0 + ^ (tfolMUHfciMifc + ll[r>,A]|U||u||y 
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where hD is a bounded operator, i. e., 

(u, A Dv)h < — ~Vq A 0 t’o + ^ (^i + A]|| & ) \\v\\\ 

Now, according to propositions 2.4, 2.5 we have L T v = 0, which is equivalent to u_ = Lv + 
(the latter L denoting the analytic boundary operator). Thus 

v^AoUo = vIA_u_ + v+A+u+ = (A+ + L t A_l) v+ > ^l u o| 2 

where the last inequality follows from the boundary conditions and the assumptions on 
L and A. Note that the analytic problem would result in exactly the same inequality. 
Hence 

(v,A Dv)h < — jl v o| 2 + g ^ ll[^’A]||/») IMIa 

Lemma 2.1 shows that 

(v,Bv) h <(\B\~ + 0(h))\\v\\l 

Consequently, 

J t \\v\\\ + |»„| 2 < m . ( 7/ ((IIIA A]|U + 2|BU + 1 + K, + 0(4))||v|g + mi) 

Integration with respect to t proves the proposition with K = max(l,2/7). □ 

Definition 3.1 A serai-discrete approximation to the initial-boundary value problem u t = 
F(x,t,d)u is said to be strictly stable, if the semi-discrete solution satisfies an energy 
estimate that is exponentially bounded by exp(a't), a' = a-\-0(h), where a is the expontial 
growth factor of the analytic estimate. 

Remark: If A (or (-,-)/i) satisfies the assumption 2.1, it follows that K\ = 0. Also, by 
lemma 2.3, \\[D, A] ||*. = p([D, A]). Eq. (16) would thus be strictly stable if p([D, A]) < 
lA'Ioo- In particular, (16) is strictly stable if A(x) = const, since this implies [D, A] = 0. 
We also point out that the proportionality constant K is completely independent of the 
discretization. In case the estimate of the boundary integral is not needed one may take 
K = 1. For variable coefficient problems we have the following result: 

Corollary 3.1 Let D and (•,•)/, be given by (3) and (4). Then (16) is strictly stable. 

Proof: 

According to the preceding remark the corollary follows if we can show that p([D, A]) < 
|A'|oo- But 

/ 0 Ai - A 0 

0.5(Ai — Ao) 0 0.5(A 2 — Ai) 

U’- M • ! 

0.5(A„_i - A„_ 3 ) 0 0.5(A I/ — A„_i) 

i A u A v _! 0 
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Assuming that A(x) is C 1 the mean value theorem gives A,- — A_, = A — j)h for some 
£ij € (ih ,jh). The corollary thus follows from the Gersgorin disc theorem. □ 


3.2 Parabolic Systems 

We consider the parabolic system 

xit — Au xx T Rxi x T Cxi F x £ (O’ 1 ) / // \ / rl \ 

u(x,0) = f(x) L 0 =( fj j L x = ( 1 ) (17) 

L o u(0,t) 4- L!U x (0,t) = 0 \ o J \ / 

where Lj, € R rflX<i , L" € R d2Xd , rfi + d 2 = d; rank(L() = rank(L^) = d 2 ; A, S, C, 

and F depend smoothly on x and t. It is assumed that the system is strongly parabolic, 
i. e., A(x, t) + A(x, t) T > 261. 

The following lemma, a proof of which can be found in [3], will be crucial when proving 
an energy estimate for the solution of (17) and its semi-discrete counterpart. 

Lemma 3.1 Let A £ R dxd be arbitrary and let L$ } L\ £ R dxd be of the form (17). The 
following conditions are equivalent: 

(i) There exists a constant c > 0 such that 

\u T Au x \ < c\u\ 2 

for all u , u x £ R d that satisfy 

L 0 u “j- L\u x — 0 

(ii) If a, be R* are vectors such that 

L[b = 0, L^a — 0 

then 

a T Ab = 0 

Assumption 3.1 Given the boundary matrices Lq, L\, the matrix A is supposed to be 
such that the second condition of lemma 3.1 holds. 

Remark: Except for Dirichlet and Neumann conditions, assumption 3.1 imposes severe 
restrictions on A. Lemma 3.1 states that the assumption above is necessary in order 
to obtain an energy estimate. The computations that follow will show how the second 
condition, which holds by assumption, implies the first. 

Before deriving the energy estimates, one more lemma is needed [3]. 
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Lemma 3.2 Suppose that assumption 3.1 holds and that A(x,t) + A(x,t) T > 261 . Then 
the d x d matrix 



is non-singular. 


As usual the boundary conditions are written as L T v — 0, where 


L t = 


Lo + d fLi 




(18) 


where doj/h are the non-zero elements of the first row of Z), which is a difference operator 
satisfying the conclusion of proposition 2.1 or 2.2. We have 


Lo 



( doo/h )I ] 



Thus, lemma 3.2 implies that L 0 + ( d 00 /h)Li is non-singular for h > 0 sufficiently small. 
From (18) it follows immediately that rank (L) = d. According to proposition 2.5 the 
corresponding projection operator is well-defined, and we obtain 


v t = P(AD 2 v + BDv + Cv + F) 

v(0) = / 


/ A(0,i) 


A = 


\ 

A(l,t) 


with similar expressions for B, C, F. 


(19) 


Proposition 3.2 Let (*,•)/! be given by (2) and suppose that D satisfies the conclusion 
of proposition 2.1. If P is defined by (10) and (18) } then the solution of (19) satisfies an 
energy estimate 

IMOIIS + j (‘ (K(r)| ! + |^(t)| 2 ) dr < e l CT '+ 0 ( A ))' (ll/HJ + j‘ ||F(r)||Jrfr) 


Proof: 

By propositions 2.5, 2.4, 2.1 we have 

(v, PAD 2 v)h = (u, AD 2 v)h = —Vq A(Dv) 0 — ( Dv , ADv)h — (u, [ D , A]Dv)h 


where we have assumed homogeneous Dirichlet conditions at the upper boundary for 
convenience. Due to proposition 2.5 it follows that 



Partition Vj = u'- + v", v'- G kerL{, v'f G (kerZ/[) x . Eq. (20) implies L^v 0 = 0 and by 
construction L\(Dv') o = 0. Hence, according to assumption 3.1 

-vlA(Dv) o = -u^(£>u") 0 

Eq. (20) can be rewritten as 

( ^ ) (Dv") o = -L 0 v o 

Since (Dv") 0 G (kerL() J_ we get 

( L ' 

T 

Li(Dv")o = -L 0 V 0 , Li = * 

K s T d2 

where {sj} is a basis in kerL{. Thus, L\ is non-singular, and one obtains 
~ v oA(Dv ) o = v^aL^LoVo < 7|u 0 | 2 , 7 = I 0 |oo 

This is exactly the same expression as one would get in the analytic case. Thus 

(u, PAD\) h < 7 |u 0 | 2 - S\\Dv\\l + ||[A A]|U||n|U||£>Hk (21) 

Furthermore, 

(v, PBDv)h < (|B|„ + O(h)) ||o|U||D»|U, (», PCv) h < (|C|„ + 0(k )) IMK (22) 

Finally, proposition 2.6 and the algebraic inequality yield 

^IMIJ+K| ! <(a' + 0(/ 1 ))||t,||J + ||F||J 

Integration with respect to time proves the proposition. □ 

Remark: All coefficients, except ||[D, A]||*, appearing in (21) and (22) are identical 
(modulus 0(/i)-terms) to those of the analytic estimate. Since the discrete Sobolev in- 
equality 2.6 introduces the same growth rate as the analytic Sobolev inequality, it follows 
that (19) is strictly stable if we have the estimate ||p, A]||/i < |A'|oo, which is true if 
A(x) = const. For variable coefficients one can prove 

Corollary 3.2 Let D and (•, )^ be given by (3) and (f). Then (19) is strictly stable if A 
is symmetric. 

□ 



Proof: 

Same as for corollary 3.1. 
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3.3 Hyperbolic-Parabolic Systems 


Consider the mixed hyperbolic-parabolic system 

U t = Au xx + B\\U X + B\ 2 V x 4" C\\ u + C\ 2 V 4- F x E (0, 1) 

Vt — Au r -|- B2\U X 4- C 21 U 4" C 22 V 4- G u E R ',t) £ R ! 


Liu*(0, 0 + L o u(0, 0 4- M o u(0, t) = 0 
v_(0,t) = 5ow+(0, t) 4- Roii(0,t ) 

u(x,0) = /(x) 
u(x,0) = <^(x) 

(23) 

As usual we assume u = v = 0 in a neighborhood of x = 1 for convenience; Lo, Ti are as 
in section 3.2, and So satisfies the hypotheses of the boundary operator in section 3.1. The 
coefficient matrices and the forcing functions of the differential equations may depend on 
x and t. 


v _ E R d i ,v+ 6 R 


d<’ 


Mo = 


Mq 


The discretized boundary conditions are written as L T w = 0, where L T E FT* x ^ +1 ^, 
d = di 4- d 2 , d' = d\ 4- d 2 is given by 

\ / d.nt \ / dn. \ \ 


r d 0 o r 

Lo 4 — Mo 

-R 0 ( / -So ) 


t*> « 

0 0 


0 

n 

0 0 


0 ... 0 


We want to show that L T has full rank. The first block of L T can be rewritten as 


/ 

(24) 


D(h) 0 
0 / 


L 0 0 \ ( ( h/doo)L ( h/doo)Moo ( hj doo)Mo\ \ 

-Ro I 0 ) + \ 0 0 -So ) 


where 


m = ( ^ 


L = 


L[ 

Li 1 




Since L is invertible, it follows that 

(l O' 
V -Ro I 


+ 


( h/doo)L (h/d 00 )Moo 
0 0 


is invertible, i. e., has full rank for h > 0 sufficiently small. The expression enclosed by 
the square brackets thus has linearly independent rows, which in turn implies that the 
first block of L T has full rank. Hence, L has full rank, and the corresponding projection 
is well-defined. 


The semi-discrete system is formulated as 

w t = P (AD*w + ADw + Cw + F) (uj\ j = 0 (25 ) 

117 ( 0 ) = \ v > J 
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Proposition 3.3 Let (-,-)h be given by (2) and suppose that D satisfies the conclusion 
of proposition 2.1. If P is defined by (10) and (24), then the solution of (25) satisfies an 
energy estimate 

IM*)lfi+ll®«lli+ £ f (KMI 2 + l»iW| 2 ) dr 

}= 0,U J0 

< Ke ia'+<Hh))l + Wi| 2 + £ (||F( T )||J + ||G(t)||J) <Jt) 

Proof: The energy method applied to (25) yields 

d 

J t \Mh = 2(u>, P(AD 2 w + A Dw + Cw + F))h = 2{w, ( AD 2 w + A Dw + Cw + F)) h 
Now 

(w,AD 2 w) h = h £ a.,(u T V T)( A ’ ° ) 1 £ d, h d u w, 

0 \ / n kj = 0 

v i v 

b 'y ' Aj —— ) ) djkdkiui — {u,AD u)h 

i,j=o n k,l = o 

where D is the difference operator of (19). The remaining terms are handled in a similar 
manner. One has 
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(i) (u>, AD 2 w)h = (u, AD 2 u) h 

(ii) ( w , kbw)h = (u, B n Du)h + (u, Bi 2 Dv)h + (v, B 2 \Du)h + {v , A Dv)h 

(iii) (w,Cw)h = ( u , C n u) h + (u, C\ 2 v)h + ( v , C 2 \u)h + (n, C^v)* 

(iv) (w, F)h = (u, F)a + (u, <?)fc 

For convenience we use the same symbol D to denote the difference operators acting on 
u and v. As far as the energy estimate is concerned, the hyperbolic-parabolic system has 
now been reduced to the previously treated hyperbolic and parabolic systems. 

Items (iii) and (iv) consist only of lower order terms, and can be estimated using 
lemma 2.2. Thus, the coefficients of the estimates are identical to the corresponding 
analytic estimate (modulo (9(/i)-terms). In item (ii) the potentially “dangerous” terms 
are those containing Dv. Using exactly the same technique as in the proof of proposition 
3.1 we get 

(v,ADv)h < — j|no| 2 + |S'jA_i?o| 00 |uo||uo| + — l-^o ^--^oloolwol 2 

+ i(tfi + ||[AA]IU)IMII 

i. e., by means of the algebraic inequality 

(v, A Dv) k < -||i>o| 2 + 7il u o| 2 + \ (Ki + l|[», A]|U) IMII 


Furthermore, 

(u, B\ 2 Dv)h < — |5i2|oo (eibol 2 + 1 l u o| 2 ) + ||[£^ #i2]IUIMUIMU — {Du, B\ 2 v)h 

Finally, in item (i) the term ( u,AD 2 u)h is treated as in the proof of proposition 3.2, the 
only difference being that 

— Uq A{Du) 0 = ulA'L^LoUo + ulAL^MoVo < 72 [e2|^o| 2 + ( e 2 1 + l) l u o| 2 ] 

We point out that the coefficients of the boundary terms in the inequalities above are 
identical to those of the analytic estimate. Choosing e\ and t 2 sufficiently small we thus 
arrive at 

J t \ Ml? + J (|u„| 2 + kl 2 ) < (o' + 0(h)) Ikii 2 + ||F||J + IIGHJ 

where we have used ||to||£ = ||u||j[ + ||v|| 2 ; in the right member we have used proposition 
2.6 and the algebraic inequality to eliminate |u 0 | 2 and ||Du|| fc . Integration proves the 
proposition with K — max(l,4/7). D 
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Remark: In case no estimate of |u 0 | 2 is needed one may take K — 1. Also, only the 
coefficients ||[Z), A]\\^ ||[Z), A]||*, \\[D, .612 ] | and K\ will be larger than their analytic 
counterparts. If either of the conditions of assumption 2.1 is met, then K\ = 0 and the 
operator norms can be replaced by the corresponding spectral radii (cf. lemma 2.3). In 
particular, if A, A, B\ 2 are constant, then (25) is strictly stable. As before, for variable 
coefficients we have 


Corollary 3.3 Let D and be given by (3) and ( f )- Then (25) is strictly stable if A 
and B u are symmetric. 

Proof: 

Same as for corollary 3.1. □ 


3.4 Strict stability 

So far we have obtained strict stability under special circumstances, such as constant 
coefficient problems or second order methods. The crux of the matter lies in estimating 
the commutator [Z), A\. Only in the previous cases were we able to prove that ||[ J D,A]||/ l < 

| Ad | oo . In fact, numerical experiments show that ||[£),A]||/ l > p{[D, A}) = A"|A/|ocn K > 1? 
for high-order methods. Typical values for ZTs corresponding to diagonal norms are 
K — 1 . 67 , K — 2 . 55 , and K = 35 . 8 , where the operator accuracy increases from three to 
five. One would still obtain K > 1 even if one considered only the interior operator. This 
indicates that the commutator should be avoided, which can be achieved if the analytic 
problem is reformulated. 

The hyperbolic system (14) can be rewritten in skew-symmetric form as 

Ut = “(A u) x + —Au x + (.B — u + F x £ (0, 1) 

u(z,0) = f(x ) 

u_(0,<) = Lu + (0,t) LeR d '* d2 

The corresponding semi-discrete system becomes 

vt = P (f DAv + \ ADv +{ B ~ ^ A ') v + F ) ( 26 ) 

u(0) = / 


Proposition 3.4 Let (•,•)/! be given by (2) and suppose that D satisfies the conclusion 
of proposition 2.1. Define P by (10) and (15). If either A or E fulfills assumption 2.1, 
then (26) is strictly stable. 
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Proof: 

The energy method implies 

4 \\v\\l = -v%Aov 0 -(Dv,Av) h + {v,ADv) h -(v,A'v) h 
at 

+2 (v,Bv) h + 2 (v,F) h 

The boundary terms are treated exactly as in the proof of proposition 3.1. Because of 
corollary 2.1 we have ( Dv,Av)h = (v,A Dv)h- Thus, by lemma 2.2, 

^IMI/1 + 2 l u °| 2 — ( l j ^ , |°° + moo + 1 + 0(h))\\v\\l + ||F||£ 

which is identical (neglecting C^(ft)-terms) to the analytic estimate. □ 

Remark: If £ is diagonal, then the 0(/i)-terms vanish identically (corollary 2.3). 

The parabolic system (17) is altered in a slightly different manner. The modified 
system reads 

u t = ( Au x ) x + (B — A')u x T Cu T F x G (0, 1) 
it(x,0) = f{x) 

L o u(0,t) + Liu x (0,t) = 0 

which is discretized as 


v t = P ( DADv + (B — A') Dv + Cv + F ) 

v(0 ) = / (Z<) 

Proposition 3.5 Let (•,•)/, be given by (2) and suppose that D satisfies the conclusion 
of proposition 2.1. If P is defined by (10) and (18), then (27) is strictly stable. 


Proof: 

Left to the reader. □ 

Finally, the mixed hyperbolic-parabolic system is reformulated as 


Ut — (Au x ) x + (B u — A'}u x + (i?i 2 v )x + C\\U -f- (C*i 2 — B[f)v + F x £ (0, 1) 
v t — -(Av) x T -Au x + Bi\u x + C 21 U + (C 22 — 2 ^') v + a 


where the initial data and the boundary conditions are identical to those of (23). In 
semi- discrete form we have 


w t = P ( DADw + DAw + BDw + (C — A')w + F ) 
ic(0) — ip 


(28) 
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where 

/ / 0 B 12 (0,<) \ 

V 0 A(0,i)/2 y 

A - 

/ 0 b 12 (m) \ 

v \ o A(i, <)/2 y 

/ Bn (0,<) -^'(0,0 0 \ 

V B 2l (0,t) A(0,<)/2 / 

5 = 

/ o A 

\ ^ B 2l {\,t) A(l,<)/2; 

Proposition 3.6 Let (vH given by (2) and suppose that D satisfies the conclusion 
of proposition 2.1. Define P by (10) and ( 24 )• If either A or E fulfills assumption 2.1 , 
then (28) is strictly stable. 

Proof: 

Left to the reader. Q 



4 Homogeneous Boundary Conditions in Two Di- 
mensions 

The results of section 3 will now be generalized to two space dimensions. If the boundary 
is smooth, the original problem can be decomposed into two problems via a partition of 
unity, one of which is a Cauchy problem. The second problem is an initial-boundary value 
problem that is periodic in one space dimension, see figure below. 



Consequently, summation by parts is needed only in one dimension, and the generalization 
of propositions 3.1, 3.2, 3.3 to two dimensions follows immediately. For details on the 
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decomposition we refer to [3]. The situation is different if the boundary is non-smooth, 
which is the case in the presence of corners. As mentioned at the end of section 2.1, it 
is not known how to extend norms of type (2) so as to obtain summation by parts in 
several space dimensions. We thus limit ourselves to diagonal norms, in which case we 
have proposition 2.3. 

All boundary conditions considered so far are local. In case of characteristic and 
Dirichlet conditions no new difficulties are presented in two dimensions, because each 
boundary point can be treated individually. Boundary conditions involving derivatives 
increase the complexity significantly. Therefore, we shall only allow normal derivatives in 
the boundary operator. This is no serious restriction from the application point of view. 
Thus, away from the corners these boundary conditions are locally one-dimensional. For 
each such boundary point we obtain a projection operator of the previous section. In 
particular, these operators commute since they affect disjoint sets of grid points. At cor- 
ners the situation is more complicated, because there are two different normal derivatives, 
which implies that the corresponding projection no longer is locally one-dimensional. 



Pi, P 2 , P c commute 


P 2 


Throughout this section we shall focus our interest on the origin, and assume that the solu- 
tions are supported only in a neighborhood of (0,0). The remaining boundary conditions 
will be accounted for by applying the projection operators corresponding to the boundary 
point in question. Since these operators commute, the resulting product is the uniquely 
defined boundary projection. The domain of definition is taken to be fi = (0, 1) x (0, 1) 
with boundary T. It will be shown later how to extend the results to curvilinear domains. 
In order to simplify the presentation all lower order terms will be omitted. 


4.1 Symmetric Hyperbolic Systems 


Consider 


2 

u t = £ A tUx , + p, *en = (o,i)x(o,i), 

t=l 

u(x,0) = f(x) x = (xi,x 2 ) 

<pi{x,t) = S{x)<pu(x,t), iff 


(29) 
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where ipj, ifjj denote the locally ingoing and outgoing characteristic variables; A,- = 
Ai(x,t), i = 1,2 are symmetric and S(x) is assumed to be “small”. It should be noted 
that 6 R dl ^, ipu € where dj(x) + d 2 (x) = d, x G T. The matrix 

2 

A(x) = ^n,(x)A,(x) (30) 

»=1 

can be diagonalized for every x 6 T; n(x) = (ni(x), 712 ( 2 )) is the outward unit normal of 
T. Hence 

A(x) = Q t (x)A(x)Q(x ), x G T (31) 

The characteristic variables are only needed at the boundary, and they are defined as 
c p{x,t ) = Q T (x)u(x, t). It will be assumed that A(x) is uniformly non-singular for x G T, 
i. e., the eigenvalues are bounded away from zero. However, the number of positive 
and negative eigenvalues may differ from one boundary point to another. The analytic 
boundary conditions can thus be expressed as 

L(x)u(x,t) = 0 L(x) = ( Qj(x) -S(x)Qjj(x) ) (32) 

Clearly, L(x) has full rank for every x E T. Strictly speaking, Z/(0,0) is not defined so 
far, because the normal n(0, 0) is not well-defined. It will soon be shown how to define 
L(0,0), and we can formally consider L(x) as being defined for every x 6 T. 

Let i = j = 0, . . . , v 2 be a grid function. Define v T = (v$ . . . uj 2 ), 

vj = [vq 3 . . . v^j). The discretized boundary conditions are written as 

LjjVj = 0, i = j = 1 , . . . , i/ 2 — 1 and j = 0, u 2 , i = (33) 

where 

Ljj = ( 0 ... 0 L(ih u jh 2 ) 0 ... o)eR J,(, ' j)x(vi+1)i 

with the non-zero element being the i th entry. At the origin we define 

L( 0, 0) = ( Qj(0, 0) -£(0, 0)^(0, 0) ) (34) 

where <5(0,0) fulfills 

2 h h 

Q T AooQ = Aoo, Aoo = ^ n,A,(0, 0), n\ — — — ,n 2 = — ; r > ^ = \/h f + h% 

i = i A h 

The motive for defining T(0,0) this way will be evident later. Furthermore, Aoo is sup- 


posed to be 

non-singular. Let 


To = 

o 

o 

t-1 

o 

6R («+i).«» i So = jr<k(i, 0 ) 

Lj = 

( ^Oj List j j 

eR («+>Hx. >j Sj= £ = 1 >/ 2 -l 

j=0,iyi 

J'l 

Ti/ 2 ~ 

( U„ 2 • * * L l/ll/2 ^ 

ER^+D^x^, = £>(;, i^) 

t=0 
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The boundary conditions may thus be expressed as 

( U \ 

L T v = 0, L = 


V 


l*2 J 


*2 


6 S = Y, S J 

j = o 


(35) 


Obviously rank(L) = s, i. e., L has full rank. Hence, the corresponding boundary projec- 
tion is well-defined, and is given by 

P = I- IT 1 L(L t £- 1 L)- 1 I t 


where 


£ = 

( <r 0 £i 

\ 

, El = 

1 ( 7 0 1 

\ 



< 7*2 Si ) 



/ 


, I G R 


dxd 


It is possible to simplify the expression for P in this case. We have 


£ L = 


But £j 1 Lj = LjHj, where 


H,= 
#i = 


If*o 

( ^M> 


( sr'io/^o \ 

Ei L^jcr^ J 

eR JjX,) j = l,. .. ,^ 2 - 1 
€ R SjXSj j = 0,1/2 


V 0- 


l-'l 


Hence 


E” 1 ! = LH, H = 


l/cr ui 

( Ho/<to 


\ 


€ R 


5X5 


Hv 2 I &v 2 ) 

Clearly, H is invertible. We therefore arrive at 

P = / - LH{L T LH)~ l L T = / - L(L t L)~'L t 

i. e., P is independent of £. 

The semi-discrete system can now be defined as 

v t =p(i: a ' d ' v + 

u(0) = / 

It will next be shown that the solution to the system above satisfies an energy estimate. 


(36) 
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Proposition 4.1 Let (-,-)h be given by (5) and suppose that D\ and D 2 satisfy the con- 
clusion of proposition 2.3. If P is defined by (10) and (35), then the solution of (36) 
satisfies an energy estimate 

IMOIfi + jf Ih’MIlr*' < *>"'* (ll/llJ + £ IIJ'MIII*-) 

where the boundary energy || • ||r is given by (v\ — v<i — v for convenience) 

IM r )llf = h 2 a 3 (Kf + kj| 2 ) + (l u *°| 2 + M 2 ) 

j = 0 »=0 


Proof: 

From propositions 2.5 and 2.4 we obtain 


— I 
dt 


Ml l = 2 A * D i V )h + 2 (*b F )h 

t=l 


From proposition 2.3 and corollary 2.1 it follows that ( v is only supported in a neighbor- 
hood of (0, 0)) 

1 ( u 

(v,AiDiv) h = -- h 2 ^2 ° 3 v1 3 A x v 0j + (v,[D 1 ,A 1 ]v) h 

1 \ j=o 

\ ( V \ 

(v, A 2 D 2 v) h = -- ( hi <rivJ 0 A 2 v i0 + (v, [ D 2 , A 2 ]v) h 

1 \ l= o / 

Thus, by lemma 2.3 we have 


d_ 

dt 


— 12 


\l < -hi °i V Jo A 222,0 -h 2 Y£ <7jV t 
»=0 j= 0 


AiVqj + Ai]) + lj \\v\\ 2 h + \\F\\ 2 h 


In the first sum the outward unit normal is n = (0,-1), and in the second n = ( — 1,0). 
Except for the origin, the boundary terms are of exactly the same form as in the one- 
dimensional case. Eqs. (30), (31) thus imply that 

~ v Jo A 2 v io = V^Ato^to < — ■~ - |</2 0 t| 2 = ^~l u oi| 2 5 i > 1 

and a similar inequality holds for the other terms. At the origin we get 

— h 2 cr qv£qA\Vqq — h\(Tov£ 0 A 2 voa = ^o^oo^oo^oo 5; —ha o-^-|<*?oo| 2 


But h > {hi + h 2 )j\/ 2. Hence 


-h 2 (T 0 v£ 0 AiV oo - hi(T 0 v£ 0 A 2 v 0 o < 


hl(T ° 2v^l V °°| 2 
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Since A(x) is uniformly non-singular it follows that 7 = inf(7oo/\/2, 7io, 7oj) > 0- Because 
of 700, the constant 7 will in general be smaller than the corresponding constant of the 
analytic energy estimate. We thus arrive at 

j t M \ i + |imi? < (e (Kia, aj) + 1) ihi; + iinu 

which proves the proposition ( K — max(l,2/7)). n 


4.2 The Heat Equation 

The analysis of a homogeneous Dirichlet condition is straightforward, even if the domain 
of definition fl is non-trivial. The problem lies in discretizing the Neumann conditions 
properly. This was clear in one space dimension. In two dimensions the occurrence of 
corners certainly complicates the analysis. To gain insight we shall begin by looking at a 
simple model problem. 

The two-dimensional heat equation reads 

lit — ^xixi T ^X2X2 X ^ (0, 1) X (0, 1) 

u n (x,t) = 0 x € r 

u(x,0) = /(x) 

where u n is the normal derivative of u. Again, we focus our attention to a neighborhood 
of (0,0). The boundary conditions are discretized as 

T~ Yl dokVhj = 0, j = 0, . . . , r d okVik = 0, i = 0, . . . , r (37) 

k= 0 k=0 

or, equivalently, 

(Div) oj = 0, j = 0, . . . , r (Z) 2 v),. o = 0, i = 0, . . . , r (38) 

where D\ and D 2 are defined by proposition 2.3. The conditions above imply that two 

boundary conditions are prescribed at the origin for the discrete problem. This approach is 

natural from the intuitive point of view, in that gradients at the origin may be interpreted 
as one-sided limits from the interior. For the time being we ignore this technicality. It 
will later be shown how it can be overcome. When deriving the projection operator it is 
convenient to cast the boundary conditions into yet a another form. Define the boundary 
operators L\j and T 2t through 

LfjV = (Div) oj = 0, j = 0, . . . , r L^v = (D 2 v) i0 = 0, i = 0,...,r (39) 

where 

Ll={ 0 ... 0 0 ••• 0 ) eR"l , ’ , ’ l|, ’«),j = 0,... 1 r 

Ll = ( ej ... 0 ... oj €R lx( " +, )(^ +,| ,i = 0 ,...,r 

\ U2 h<2 / 
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Here {e,} is the canonical basis in R l ' l+1 . The boundary conditions can thus be written 
in standard form 

L T v = 0 L=(L 10 ... L\ r L 20 ... T 2r ) € R 2(r+1)x(l,1+1)( ^ +1) (40) 

We know that the corresponding projection operator is well-defined iff rank(T) = 2(r-f 1). 

Lemma 4.1 The columns of L (40) are linearly dependent. In particular, rank(L) < 
2r + 1 . 

Proof: 

To investigate linear dependence we study 

r r 

E a jL\j + E PjL2j = 0 

3=0 j = 0 


which is equivalent to 


E (ctjh 2 d 0 k + 0khid oj ) e k = 0, j = 0, . . . , r 


k = 0 


Since {ejJ is an ortho-normal system it follows that 

Otj h2(%ok “h $kh i ^Oj J 5 ^ 0, . . . , r 

which obviously has the non-trivial solution 


OLj C^Ojr ft j (Iq j ? J . 5 


h x 


and the lemma is proved. 


□ 


As a consequence of lemma 4.1 the projection formulation breaks down. If, however, 
we change the boundary condition at the origin to 


L lx v = ((' “ X)ito + xilo) » = 0 0 < X < 1 


(41) 


and leave the boundary conditions at the remaining points unchanged we get a well-defined 
projection operator, since 


L=(L u ■■■ L u L 0x Ln ... L„ ) 6 


(42) 


has full rank. 


Lemma 4.2 The columns of L (42) are linearly independent. In particular, rank(L) = 
2r + 1. 
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Proof: 

Again we study 

r r 

T, <*jL\j + 7 L 0x + y PjL 2] = 0 

i=i i= i 

which is the same as 


dm £ fok + 1 f(l - x)r £ + x^e„) = 0 

‘2 fc = l V h l k = 0 ft 2 / 

r 1 r 


77 yfok + OjT" yd ok e k + 7X7^0 = 0 j = 


2 Jb=l 


iiS 


The first component of the first equation yields 7(L 2 (1 — x) + ^iX)^oo = 0. Since doo ^ 0 
for any operator satisfying proposition 2.3, and since L, > 0, 0 < x ^ 1) necessarily 
7 = 0. From the remaining components of the first equation we then obtain j3j — 0, 
j = l,...,r, which in turn implies aj = 0, j = l,...,r. The columns of L are thus 
linearly independent, i. e., L has full rank. □ 

Before proceeding with the energy estimate, one more technical lemma is needed. Let 
Lo x , and Lq X2 be defined by (41), and let 

L=(Ln ... L u L 0xi U X2 L 21 ... L 2r ) e R 2(r+1)x(l/1+1) ^ +1) (43) 


Lemma 4.3 The columns of L (43) are linearly dependent. In particular, rank(L) < 
2r + 1. 


Proof: 

Consider 

T T 

y oijL\j + 71 L 0xi + 7 2 Lo X2 + y fl]L 2j = 0 

j=i i=i 

Obviously the lemma is true for = \ 2 - In the following we thus assume xi i 1 X2- The 
equation above can be rewritten as 


T 


r 


y a jLij + y 0jL 2 j 

j = o j=o 


= 0 


where 



According to lemma 4.1, eq. (44) has the non-trivial solution 


atj — doj 


ft = ~doj, 

fix 


J = 0, . - . , r 


(44) 
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whence 


71 = 

72 = 


X2 + 7— (1 — X2) ) /(x 2 — Xl) 
Xi + 77(1 — Xi) ) /(x 2 — Xi) 

'll / 




solves the original equation. The lemma is proved. 


0j = —rd oj, j = l,...,r 

'ii 


□ 


Proposition 4.2 Let P 6e given by proposition 2.5, where L is defined by (42). Then 
Lj 0 P = L^qP = 0. 


Proof: 

Clearly, L T P = 0. Furthermore, L w , L 2 o = L 0x for x = 0,1, respectively. But then, by 
lemma 4.3, 

L w = Lat\ L 2 o — Lot2 

for some vectors 01,02 G R 2r+1 . This proves the proposition. □ 

Remark: Suppose that v is a vector such that v = Pv , where P is as in the previous 
proposition. Then L w v = L 20 v = 0, i. e., (Piu)oo = {D 2 v) 0 q = 0. In other words, by 
requiring that the boundary condition at the origin hold for a specific convex combination 
we actually get the stronger result (Z)it>)oo = {D 2 v)qq = 0- Thus, we need not overspecify 
at the corners, cf. eq. (37). In the appendix we give a direct proof that Lj 0 P = 0 for L$ x 
with x = 0.5. 

The semi-discrete heat equation is given by 

v, = P(D\ + Dl)v 

»(o) = f (4a) 

Proposition 4.3 Let be given by (5) and suppose that D\ and D 2 satisfy the con- 

clusion of proposition 2.3. If P is defined by (10) and (42), then the solution of (45) 
satisfies an energy estimate 

IKOIU <11/11* 


Proof: 

The energy method gives 

j t \\v\\l=2(v,Djv) k +2(v,Dlv) h 

By proposition 2.3 (v is supported only in a neighborhood of the origin), 

r 

(v,D*v) h = -h 2 '^(JjVQi(D l v)o ] - ||T>iHIa 
j=o 

According to propositions 2.4, 2.5 and 4.2 we have 

(Thu)oj = LjjV = 0, j = 0, ... ,r 

The remaining term (u, D\v)h is treated similarly, and the proposition follows. 


□ 
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4.3 Parabolic Systems 

Consider 


Ut — ^ ] AjjU x iX) F , 

«,i=i 

u(x,0) = /(x), x = (xj , x 2 ) 

L 0 (z)u(x, t) + Li(x)u n (x, f) = 0, xGT 


x € ft = (0, 1) x (0, 1), w6R j 


(46) 


The assumptions on L 0 ,Li in (17) are supposed to hold pointwise for each x € T. Fur- 
thermore, we require that assumption 3.1 with A = An be valid on x, = 0, i = 1,2. In 
particular, the conclusion of lemma 3.2 holds for each boundary point. It will be assumed 
that (46) is strongly parabolic, i. e., there are vectors itj(x,<) € R^, i = 1,2, such that 

2 2 

J2 Ui(x,t) T Aij(x,t)uj(x,t) > 26^|u,-(x,t)| 2 

t.j=i t=i 

for all x 6 0, t > 0. If the matrices Ajj ^ 0 , j / j, then the assumptions must be 
strengthened. The energy method applied to one of the cross terms yields (u is supported 
only at the origin, An — const for simplicity, f 1 is the unit square) 

(u, A l2 u XlX2 ) = - / u T A l2 u X2 dx 2 - {u Xl , A 12 u Xl ) 

Jx 1=0 

In general we cannot get an estimate of u X2 (0, x 2 , t) in the boundary integral. It is therefore 
natural to require 


Assumption 4.1 Ajj = Aij, i ^ j . 
Remark: Neglecting scaling factors we have 

A12 = A 2 i = — 



0 

0 

0 

°\ 

1 

0 

0 

1 

0 

p 

0 

1 

0 

0 


0 

0 

0 

0 / 


for the Navier-Stokes equations ( p denotes the density). Clearly, assumption 4.1 is fulfilled. 
If assumption 4.1 holds one can integrate by parts once more to obtain 

(u, A\ 2 u XlX2 ) = — it y4j2(0,0,t)u (^xj , Ai 2 u X2 ) 

In two dimensions we cannot eliminate the boundary terms by means of Sobolev inequal- 
ities, since they would involve T 2 -norms of u XiXl and so forth. This motivates 
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Assumption 4.2 Let u(x,t) satisfy 


Z>o(0,0)'u + Li(0, 0)u n — 0 

at the origin . Then 

u T Aij( 0, 0, t)u = 0, i ^ j 


Remark: This assumption ensures an energy estimate for the continuous problem in 
case of a non-smooth boundary, and couples the cross terms of the differential operator 
to the boundary conditions at the origin. In case of the Navier-Stokes equations one has 
zero velocity at the origin. Hence, the state vector becomes u T = ( p 0 0 p ), which 
implies assumption 4.2. 

The discrete boundary conditions are formulated as (D\ and D 2 are defined by propo- 
sition 2.3) 


Ljjv = L o (0,jh 2 )v O j + Li(0,jh 2 ) (Div) oj = j = 0, . . . , r 
L\ { v = Z/ O (^i,0)fio + Li(ihi,0) ( D 2 v) i0 = 0, i = 0, . . . ,r 

where 

Ljj = ^ 0 ... 0 Lo(0,jh 2 )eQ Li(0,jh 2 )—^2^ok^J 0 ••• 0 


(47) 


Ll 


L 0 (ihi ,0) + L\ (ih\ , 0) 


MM 


M(ih„ 0)^ef 
h2 


and ej = ( 0 ... 0 I 0 ... 0 ) € R dx ( ,/ > +1 )‘ i . The boundary conditions can be 
expressed in the usual form as 


L t v = 0 L=(l n ... L\ r L 0x L 2l ... i 2r ) 6 R (2r+1)Jx(l, ‘ +Mrf (48) 


where 

Lq x = (1 — x)-^io + X^2o 0<X<1 


Lemma 4.4 The columns of L (f8) are linearly independent for sufficiently small step 
lengths hi and h 2 . In particular, rank(L) = (2 r + l)d. 


Proof: 

Imitating the proof of lemma 4.2 gives 


7 


Lq(0, 0) + doo 




= 0 
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By lemma 3.2 the expression inside the brackets is non-singular for h\,h 2 sufficiently 
small. Hence, 7 = 0, which in turn implies a 3 = (i 2 = 0, j = 1 , . . . , r. Since the columns 
of each block L\j, L 2 j and L 0x are linearly independent, the lemma follows. □ 

The semi-discrete parabolic system reads 

v t = P ( Yj AjDiDjV -(- F 
\«\j= i 
u(0) = / 

where P is defined by proposition 2.5 and by (48). Unfortunately, assumption 4.2 is not 
sufficient for the semi-discrete problem. We need 

Assumption 4.3 Let v satisfy 

Lo{0,0)voo + Li(0,Q)((1-x)(D 1 v) oo + x {D 2 v)oo) = 0 0 < x < 1 

at the origin. Then 

(i) 

u ooA?'(°> 0 u oo = 0, i ^ j 



(ii) 


UocAi(0,0,f) = uJo^22(0, 0, t ) 


Remark: The first requirement is identical to that of assumption 4.2. The second, 
however, appears only in the discrete case. We note that assumption 4.3 holds for the 
Navier-Stokes equations, since A\\ and A 22 are given by 


/ 0 0 0 0 \ 

\_ 0 cj 0 0 

P 0 0 10 

\ -C 2 pjp 0 0 c 2 / 


a 22 


0 0 0 
\_ 0 10 

p 0 0 ci 

\ -C2 p/p o 0 


0 \ 
0 
0 

C2 / 


Hence, v% 0 A u = vlf 0 A 22 = c 2 ( -p 2 / p 2 0 0 p/p ). 


Proposition 4.4 Let be given by (5) and suppose that D i and D 2 satisfy the con- 

clusion of proposition 2.3. If P is defined by (10) and (48), and if assumptions 4-1 and 
4-3 hold, then the solution of (49) satisfies an energy estimate 

IM0II* + £ I |®(r)||frfr < (||/||J + £ ||F(r)||^r) 
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Proof: 

The energy method gives 

d 2 / r r 

-d Ml < -Z^lh^VkVokAijiDjvhk + hJ^VkvloAijiDrfko 
01 j = 1 V *=0 fc=0 

-2<s£ MAt-llJ + (Xo + CW) IMI1 + Ilf IIS 

2 = 1 

where K 0 depends on \\[D{, A,-,*] Ha, * = 1,2 and p([Di, Aij]), i ^ j. The first cross term 
can be written as (v has compact support). 

— k 2 (T kVQ k A 12 {D 2 v)ok = ~h 2 ^ ( T” ^2 dklV ° l ) = — (^0, Ai 2 D 2 V 0 )h 2 

fc= 0 *=0 V^ 2 /=0 / 

where nj = (uj) . . . and where Z) 2 satisfies (1) with respect to the one-dimensional 
scalar product (*,*)a 2 - Hence, 

— (v 0 ^ Ai 2 D 2 vo)h 2 = 2^00^12(^1 0? 0 v oo + [D2, Ai 2 ]vo)h 2 

By assumption 4.3 the boundary terms vanish. The remaining cross term is treated in a 
similar vein. 

Next, we take care of the boundary terms corresponding to the pure second differences. 
Only the origin needs to be analyzed, since the other boundary points are treated exactly 
as in the proof of proposition 3.2. At the origin we get 

— h 2 aoVQ 0 Au{Div)oo — ^1^0^00^22(^2^)00 — 

— {h\ + h 2 )cr 0 VQ Q ((1 — x)An(Div)oo + X^22(T*2^)oo) ? X = ^ ^ 

and, by assumption 4.3, 

-h 2 <T 0 Vo 0 A u (Div)oo - hicr 0 Vo 0 A 22 (D 2 v) 0 o = 

— {hi + h 2 )(ToVQ 0 Au ((1 — x){D\v)oo + x(^2^)oo) 

But v = Pv implies L 0x v = 0, i. e., by (47) 

Lo(0, 0)uoo + li(0,0)((l — x)(^i u )oo + x(^2^)oo) — 0 

In particular, Ljfv 00 = 0. Partition v = v' + v" where v'- € ker(l{), v"- € ker(L()- L . 
Assumption 3.1 then gives 

— h2(ToVQ Q Au(Div)oQ — h\(ToVQ 0 A22{D2v)oo = 

— (hi + h2)cr 0 VQ 0 A u ((1 — x)(Au")oo + x(^2w")oo) 
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By construction 


Z/o(0, 0)uoo + ^i(0, 0) ((1 — xX^i'Ooo + x{D2 v ")oo) — 0 

which can be solved in exactly the same way as the corresponding equation in the proof 
of proposition 3.2. Hence, 

— h 2 cT 0 VQ 0 A u (Div)oo — h\< 7 QVQ 0 A 22 {D 2 v)oQ = h, 2 croVQ 0 AiiL l 1 LoVoo + ^ 1 ^ 0 ^ 00 ^ 22 -^ 1 1 ^ 0^00 
where we again have invoked assumption 4.3. We thus arrive at 

■jrlMlj’ + |M|p < ( 2 \AnL\ 1 ^oh.oo + p([^ 2 , ^12]) + l) ^2 ffk\vok \ 2 

dt k = 0 

r 

+ ^21^422-^1 Uo|i,oo + p([D\, ^ 21 ]) + l) hi y] ar k \v k0 \ 2 

k= 0 

- 2«f:iiD j »iu+(/fo+ow)iHij+ira 

i=l 

where |u| li00 = sup(|v*o|) and |u| 2 ,oo = sup(|u 0 jt|)- Replacing p([Di,Aji]) by |A' -| l>00 , i / j, 
one obtains the coefficients of the boundary terms of the analytic energy estimate. They 
are thus identical if the coefficient matrices are constant or if we use the standard second 
order method. Finally, the boundary terms of the right hand are eliminated by applying 
the one-dimensional Sobolev inequality 2.6 in the xj- and x 2 -directions, respectively. This 
proves the proposition. □ 

Remark: It is clear from the proof that (49) is strictly stable if the coefficient matrices 
are constant, or if A? = An, i t = 1,2 and the second-order method (3) is used. 


4.4 Hyperbolic-Parabolic Systems 

In this section we merely formulate the problem and state the main result. The reader is 
asked to fill in the details. We consider the mixed hyperbolic-parabolic system 

2 2 

^ ^ AijU XiXj T C\ iV x> T F r £ 

i,J=l >=i 

2 2 

Vt — ) ^ B{V Xt -)- ) " C% iU x , T G ^ 6 

i=i t=i 

Li(x)u n (x, t) + L 0 (x)u(x, t ) + M(x)v(x, t ) = 0 x 6 
< pi(x,t) = S(x)ipu(x,t) + R(x)u(x,t) x £ 

u(x,0) = f{x) 
u(x,0) = <f>(x) 

(50) 


n 


R dl ,u e R d2 

r 

r 


M(x) = 


M'(x) 

0 
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The structural hypotheses on the parabolic and hyperbolic parts of the principal operator 
are identical to those in sections 4.3 and 4.1. This remark also pertains to the boundary 
conditions. In particular, the characteristic variables ipi and 9 ?// are defined as in section 
4.1. The principal operator is described by the A,j’s and Bf s; the lower order coupling is 
determined by the C,j’s. Similarly, in the boundary conditions the coupling is expressed 
by M and R. 

Define w T = ( u T v T ) and let P be the projection corresponding to the boundary 
conditions of (50). The semi-discrete system is then defined as 


w t — P (£ AijDiDjW + B t D,w + J2 C t DiW + F 

1 1=1 1=1 


ru( 0 ) = 


where 


( ( ^(0,0,0 0 

l 0 0 


A{j — 


Bi = 


\ 


( ( 0 0 

lo Bi (0,0,0 


\ V 

( ( 0 Cm (0,0,0 

l Cai(0,0,i) 0 


A, ,(1,1,0 0\ 

0 0 ) ) 

\ 


0 0 \ 

0 Bi (1,1,0 j ) 


Ci = 


L/,j 


0 Cm (1,1,0 

C 2l (1,1,0 0 


(51) 


The forcing function F and the initial data V’ 


are defined in a similar fashion. 


Proposition 4.5 Let (•,•)/* be given by (5) and suppose that D satisfies the conclusion 
of proposition 2.1. Then the solution of (51) satisfies an energy estimate 


ll“Wlli + ll»MIK + f (ll u ( T )llr + ll u ( T )llr) d T 

< (ll/Uj + W|| 2 , + f (||F(r)||J+ ||G(r)|| 2 ) dr) 
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5 General Domains and Strict Stability 


Nothing has been said about strict stability in two dimensions thus far. The purpose of 
strict stability is to ensure the same growth rate of the discrete and analytic solutions. 
If the analytic problem is defined on a curvilinear domain fi, then there must exist a 
diffeomorphism £ = £(x) of fl onto the unit square (0, 1) x (0, 1) in order for the finite 
difference method to be well-defined. Consequently, a constant coefficient problem in the 
original domain may be transformed to a variable coefficient problem on the unit square, 
which may account for a non-physical growth in the discrete estimate. 


Let £ = £(x) be a diffeomorphism of f! onto / = (0, 1) x (0, 1). The following identities 
are readily established 


dx\ 


dx 2 

— j-;* 

d(i 

dx 2 

d(i 

dx\ 

dxi 

7 

\ 

II 

dx 2 


<96 

OX2 

56 

dx\ 


J = det 



which in turn implies 

t (■'"%), = 0 

«=i 


(52) 


(53) 


where d denotes the two-dimensional gradient operator. We require that £(x) be uniformly 
non-singular, i. e., there exists a constant 6 > 0 such that J~ l > 8 on f), For later use we 
record the normal and tangential derivatives u n , and u Ti at the boundaries corresponding 
to 6 = 0, i = 1, 2: 


u r , = (-l)'%/|i(,| 

un. = -(«£ i ■ + at. ■ dt,ui,)l\at,\ 


*=1,2 j 4- i 


(54) 


where the boundary T of the domain f l has been parametrized in the positive direction. 
The analytic scalar product obeys 

(u, v) = / u T (x)v(x)dr = f u T (x(^))v(x(())J~ 1 d^ 

Jn Ji 

which suggests the following semi-discrete scalar product 


(u,v) h = (u,J 1 v) h = (J 1 u,v ) l l 


(55) 


where the last equality follows since J 1 and E are diagonal. Thus, each grid point is scaled 
with the cell volume. Similarly, the analytic boundary integrals can be parametrized as 


f u T (x)v(x)ds = f u T (x(£ 1 ,0))u(x(£ 1 ,0))|x£,(£ 1 ,0)|d£ 1 
Jr, Jo 
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Hence, it is natural to define the boundary scalar product as 

V V 

{u, v)r = J2 a J ( s 0 j u 0 ] v 0 j + a ' ( s iouJo v io + SivuTvivJ (56) 

j=0 i=0 

where the arc lengths are defined as 

$oj = \x^(0,jh 2 )\h 2 s t0 = |* {l (tfc,,0)|Ai h l = A( u h 2 = A£ 2 
with similar definitions for s u j and 

In order to prove stability we must have Pv = v. Since v will be the solution of 
equations like (60), proposition 2.4 implies that PJ~ l v = J~ l v. Therefore, it is natural 
to require 

PJ - 1 = r l P (57) 

For a general P this identity expresses a compatibility condition between the analytic 
boundary conditions and the mapping £(:c). Let P be given by (10) and (42). Then (57) 
certainly holds if 

Jij s J(x(ih u jh 2 )) = Jioy j = 0,...,r 

Jii = J(x(ih u jh 2 )) = J 0j , i = 0, . . . ,r 1 j 

which states that the mapping £(x) is locally isochoric in the rr^-direction at the boundary, 
where d/d£i is the non-tangential derivative. In case of characteristic boundary conditions 
and Dirichlet conditions we have r = 0, and (58) is trivially satisfied. For general bound- 
ary conditions, however, (58) couples the boundary operator to the grid transformation 
(cf. assumption 4.3, which links the differential operator to the boundary operator). 


5.1 Symmetric Hyperbolic Systems 


Using (53) we recast (29) into a form that eliminates the need for the commutator in the 
semi-discrete case. 



\ E + ■r'fl.Uf.) - ij-'divMu + J-'F 


(59) 


where 


Bi = dii • a = 


dji_ 

dx\ 


A\ + 



div(yl) = d • A = 


dAi dM 

dx\ dx 2 


The boundary conditions are as in (29), except at the origin, where we require that 
the characteristic boundary conditions ^>/(0,t) = 5(0)</>//(0, t) be satisfied for <^(0, t) = 
Qf (0)u(0, t), i — 1,2, where Qj (O)y^(O, t)Qi(0) are diagonal, which means that two 
boundary conditions are prescribed at the corner. This situation occurs for the Euler 
equations at corners, where it is natural to require that both normal components of the 
velocity field be zero. Furthermore, this assumption simplifies the computations that 
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follow. The projection operator is still well-defined. It should be pointed out that the 
boundary condition at the origin is only used for the semi-discrete system. 

Eq. (59) is discretized as 

(J~ 1 v) ( = P Q £ + J- x BiDiv) - \j~ l Cv + J- 1 ^) (60) 

with 

div(,4)(x(0, 0), t) 

div(i4)(x(l, 1), i) 

Proposition 5.1 The approximation (60) is strictly stable . 

Proof: 

The energy method yields (using P(J~ l v) — c/ _1 u, PJ -1 = J” 1 P => Pv = v) 

^(v, v) h = £ ((v, DiJ- x Biv)h + (v, J- x BiDiv) h ) - («, J-'Cv) h + 2(v, J~ l F) h 
at i = 1 

But (using Pi(0, jh 2 ) instead of Pi(x(0, j/ 12 ), 0 and so forth to make the notation less 
cumbersome) 

(u, D\J~ X B\v)h = -h2^2 <T j v ojJ~ 1 (^Jh 2 )B 1 (OJh 2 )v 03 - (D x v, J~ x B x v) h 

j = 0 

Since diagonal scalar products are used, assumption 2.1 holds a fortiori. Hence, (Bf = 
Bi) 

( DiV,J 1 B x v) h = (B X J 1 D x v,v) h = (J 1 BiDiv,v) h 
where the last equality follows since B\ and J~ x commute. Thus, 

V 

(v,D x J~ x B x v) h = ~h 2 '52 cr j v o : J- 1 (Q,3 h 2)Bi(0,jh2)voj - {J~ l B x D x v, v) h 

3 = 0 

with a similar relation for (u, D 2 J~ 1 B 2 v) k . We thus arrive at 

d v 

■J {v,v) h < -/i 2 ^cr i uJj _1 ( 0 ,i/i 2 )Bi( 0 ,j/i 2 )voj 

at j=o 

-hi^2 (TivJ 0 J~ l (ih x , 0) B 2 (ih x , 0)u io 

1=0 

+(|div(i4)|oo + !)(«,«)* + {F,F) h 
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By means of ( 52 ) it follows that 


= — -A,-— '-A, 


Apparently, x ^ 2 is a tangent vector of the curve x(0,^2)- Hence 

/ dx 2 dx-j \ 

V <96 #6 ) 

is an outward normal to x(0,^2) € f. The unit normal is then defined as 


ni n 2 


)-(-S si' 1 ' 


Using the definition of the arc length we then obtain 

V V 

-h 2 '}2<7 J vl J J~' l {0,jh 2 )B l (Q,ih 2 )v Oj = ^2 (TjSojVqj (nxA x + n 2 A 2 ) v 0j 

3=0 j = 0 

The boundary conditions are satisfied, whence 

vL ( nx Ax + n 2 A 2 )v 0j < ~lj\voj\ 2 


Letting 7 = inf (7.7) > 0 implies 


—h 2 ^ &j v ojJ '(OJhJBxiOJhJvo, < -7^2(TjSoj\voj\ 2 

3 = 0 j = 0 

We have thus established 

+ l{ v i v )r < (|div(A)|oo + l)(v,v) h + {F,F) h 

This is exactly the same estimate that one would get in the analytic case, and the propo- 
sition has been proved. □ 


5.2 The Heat Equation 

The heat equation in self-adjoint form reads 
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At the boundary the normal derivative is set to zero. Define 


1 i« 0 ' 0 » 


M tJ = 


96 

dx t 


(*(1,1)) 


/ 


and 


A = £ M tJ D } 

3 = 1 


Clearly, Z), is a consistent approximation of 9/9a;,. Let 


-Loni 


-|96|Tio- 

“ 1 961620 — 



(62) 


be approximations of the normal derivatives at the origin. The boundary projection P is 
defined by (10) and (42) with L 0x replaced by either of L 0 n , and L 0n2 . 

It should be noted that P may no longer be unconditionally well-defined. Arguing 
exactly as in the proof of lemma 4.2 one obtains (using Lq x = Lo ni ) 


-(|96l% + 96- 96 * 1)7 = 0 

if 96-^6 < 0, i. e., at acute corners there is a possibility of a non-zero 7 if 


h 2 


96 • 96 
|96l 2 


6 


(63) 


Hence, at acute corners we assume that h\ and h ,2 be such that (63) does not hold. 
Furthermore, lemma 4.3 is valid for L 0xi = Zo n , , Z/o X2 = 6)n 2 - This is obvious if 96 -96 = 
±|96||96|, because then L 0ni = ±L o n2 - Otherwise, we obtain (44) where 


1961 
96 • 96 
V 1961 


96 ■ 96 \ 

1961 ( 7. \ 

1961 



which has a unique solution iff |96 • 961 < |96||96l- 

The semi-discrete heat equation is now defined as 



Proposition 5.2 Assume that the mapping £(fi) = / is locally isochoric at the boundary 
in the sense of (58), and that the grid is orthogonal at the boundaries except at the corners. 
Then (64) is strictly stable. 


Proof: 

Since the transformation is locally isochoric at the boundary we get Pv = v. Thus, the 
energy method implies 

d 2 

— (v,v) h = 2 £(v,D fc J 1 M ik D i v) h 

i,fc= 1 

Summation by parts yields (u is assumed to have compact support) 


dt 


(v,v) h = -2 £> 1 M il Div)oi - 2£ hi £ <T t vi 0 {J 1 M i2 D t v) 


10 


*=1 /=0 


i=l /=0 


-2 Y^{ M ikD k v,J l Div) h 

J,fc=l 

Obviously, 

2 2 
£ (M ik D k v,J- x Div) h = Y,(DiV,biv) h 

i,k= 1 1 = 1 

Next we turn our attention to the boundary terms. We have 

E(J = Vl^iloi f 1^6 |o/(^i^)o/ + ^ ^ (^ 2 ^) 0 / 


1=1 


\d£i\oi 


The parenthetical expression is recognized as a discretization of the normal derivative 
(cf. (54)). The other boundary is treated analogously. At £1 = 0 we thus define a “normal 
difference” operator D ni through 

(D ni v)oi = — (\d£i\oi{Div)oi + - 5 “(^ 2 ^) 0 / 
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with a similar definition of D U2 at = 0. From (52) it follows that J 0 / 1 |^i| = 
Hence, 

d - 2 

^ {v } v) k = 2 (u, P n u)r - 

Using v = Pv and the orthogonality assumptions it follows that 

(D ni v)oi = = 0 (. D n 2 v)io = -|^2|/o^2/ v = 0 / > 0 


At the origin we have 


(D nt v)oo = Ll n V = 0 (Z> n2 u)oo = = 0 
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where Lq^v vanishes because of the construction of P; Lq U7 v disappears since we have 
shown that Lo n2 belongs to the column space of P (cf. proposition 4.2). Hence the 
boundary sum is identically zero, which proves the proposition. □ 

Remark: It would still be possible to prove strict stability, even if the grid were not 
orthogonal at the boundary. To compensate for the loss of orthogonality it is necessary 
to require that the grid be globally isochoric in a neighborhood of the boundary T. 


5.3 General Parabolic Systems 


When considering parabolic equations in general, tangential derivatives may appear in the 
boundary integrals, potentially calling for integration by parts once more. The occurrence 
of tangential derivatives depends on the coefficients of the original equation, the geometry, 
and the presence of mixed derivatives. These criteria are not independent of one another. 
The following simple example will illustrate this interdependence. Consider the parabolic 
model equation 

lit — U xjxi T W x\x 2 T ^X2^2 x G (b5) 

where Q is diffeomorphic to the unit square; the boundary conditions are yet to be spec- 
ified. The energy method gives (the cross term is integrated with respect to xi) 


7 ||^|| — 2 / (tzii n Tt\ xixi X 2 )ds 2 f T -b zx^tZxjjdx 

at Jr Ju 

The normal and tangential derivatives are defined as 


d 

d 

d 


d 

d 

d 

dn 

II 

+ n 2 

OX 2 

4=^ 

dx\ 

= 

+ Tl d~r 

d 

d 

d 


d 

d 

d 

dr 

t-h 

N 

II 

+ r 2 0 

UX2 


dx 2 

= n ^ 

+ T 2 — 
OT 


where n is the outward unit normal as usual; the unit tangential r is chosen corresponding 
to a positive orientation of T. Thus 


Ti = -n 2 

r 2 = rii 


( 66 ) 


If, on the other hand, the cross term is integrated with respect to x 2 , we obtain 

-IMI 2 = 2 [ (uu n + n 2 uu Xl )ds - 2 / (u Xl u Xl + u Xl u X2 + u X2 u X2 )dx 
Jr Jo 


d_ 

Jt' 


We must show that 


J niuu x2 ds = J n 2 uu Xl ds 


(67) 
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in order for the energy method to be well-defined. Using the definitions above gives 

J n^uu X2 ds = J (nin 2 uu n + riiT 2 uu T )ds 

J n 2 uu Xl ds = J (riin 2 uu n + n 2 r-[UU r )ds 
Clearly, (67) will follow iff 

J n\T 2 uu r ds = J n 2 Tiuu r ds 

From (66) and n\ + n\ = 1 it follows immediately that 

J niT 2 uu r ds = J n 2 riuu r ds — J uu r ds 

Note that the second integral of the right hand side would vanish identically if T were 
smooth. To simplify the analysis it will be supposed that u is supported only in a neigh- 
borhood of the lower left corner. Hence, it will be sufficient to consider the boundary 
portions IT and T 2 corresponding to = 0 and £2 = 0. Parametrizing T in the positive 
direction gives (cf. (54)) 


J^uu T ds = ^j^ -{u 2 )t 2 d( 2 + ^ j Q (u 2 ) f i<J£ 


1 r\ 


Letting £2 —£2 in the first integral of the right member gives (d/d£ 2 

dh -> dl 2 ) 

J T uu T ds = - T (u 2 )i 2 d£ 2 + - jf (ix 2 ) 6 = 0 


-dfdfa 


and (67) follows. The energy method is thus well-defined, and we have 
d_ 

dt" " Jr 
,2 


«-il 2 + ltoJI 2 


IM| 2 < % J {{ 1 + riin 2 )uu n + n\uu T )ds — (||i 

The quantity n 2 is discontinuous at the corners. Define the jump discontinuity 

\n\\(x) = n\ R (x) - n\ L (x) 

where n 2 fi (a:) and n\ L (x) are the left and right limits of n 2 at x (according to the positive 
orientation of T). Straightforward computations show that 

jf n\uu T ds = i ^[n 2 ](a; C| )u 2 (a: ct -, t) - ^ J^{n\) T u 2 ds 

where x c i, i — 1, . . . ,4 are the corner points. Thus, 


dt 


u 


II 2 < ]£[*?](*« )« 2 + J T (2(1 + n\n 2 )uu n - (n 2 ) T u 2 ) ds - (||u Xl || 2 + ||u X2 || 2 ) 
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From this inequality it is obvious that giving Dirichlet data at the corners and Neumann 
data at the remaining boundaries would yield an energy estimate. In fact, we could even 
allow inhomogeneous Dirichlet data at the corners and still obtain an energy estimate in 
terms of the data. The effect of the corners disappears iff [uj](x) = 0, which happens iff 


(i) n xL {x) = n lR (x) 

(ii) n lL (x) = -n lR (x) 


The first case implies that the normal is continuous, i. e., x is not a corner point. The 
second case is more interesting, since the normal is discontinuous, but the effect on the 
energy estimate disappears. This illustrates how the geometry can interact with the cross 
terms. The simplest example is obtained by solving (65) on D being the square with (1,0), 
(0,1), ( — 1,0), and (0,-1) as vertices. Evidently, the second case holds at the corners, 
and no corner values should appear in the energy estimate. This can also be seen by a 
change of coordinates: 

6 = --L I| + -L* 2 

Eq. (65) is then transformed into 


u t ~ 2 U *i<i + 


, , 1 1 X / 1 1 X 

v^’v^ )X( v / 2’V2 ) 


The cross term has vanished; instead the equation has become anisotropic. Working in 
this coordinate system it is apparent that no tangential derivatives — and consequently 
no point values — will appear when deriving the energy estimate. 

To solve (65) by means of finite difference methods it is rewritten in self-adjoint form: 


= £ (•/-’ (l + - W ) • H + D- 1 )' 

k= 1 k±l 


t k 


(68) 


where n W = — This equation is discretized in space the usual way. The cross 

terms 

(-1)' (n\ k) n[ k) u (l )^ 

must be integrated twice to eliminate the tangential derivatives. In the semi-discrete case 
this amounts to performing summation by parts twice, the second of which will require the 
introduction of a commutator, thereby obliterating strict stability (except for the second 
order accurate difference operator). To restore strict stability it would be tempting to 
reformulate the critical terms in skew-symmetric form: 


= 5 - 5 W ) 6 » 
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Doing so, however, would introduce lower order energy terms {,}/*, whose presence would 
destroy strict stability. The simplest way to resolve this ambiguity is to assume homoge- 
neous Dirichlet data, in which case the boundary terms vanish identically, and (68) would 
be the preferred choice. The choice of homogeneous Dirichlet data to eliminate the influ- 
ence of the mixed derivatives arises naturally when solving the Navier- Stokes equations, 
since at solid boundaries we have zero velocity, and since the cross terms involve only the 
velocity components. 


We now turn to general parabolic systems subject to homogeneous Dirichlet conditions. 
For general domains D eq. (46) is written as 


(j = ^2 (j 1 ^-A i jU x )\ 1 div{Aj)u X} + J 1 F 

* V ax ‘ / i k i=i 


u(x,0) = f(x) 
u(x , t) = 0 


(69) 


x e r 


where div(A,) = (Aij) Xl + ( ^2j)x 2 • Define Cj = div(v4j). The semi-discrete system is 
then given by 


( 2 2 

£ D k J~ l MikAijDjV - £ J-'CjD, v + J~'F 

i,j,k=l j=l 

The projection operator P represents the homogeneous Dirichlet conditions. 


(70) 


Proposition 5.3 The approximation (70) is strictly stable. 


Proof: 

Left to the reader □ 


6 Inhomogeneous Boundary Conditions 


The principle for handling inhomogeneous boundary data is best illustrated by means of 
a simple example. Consider the one-dimensional advection equation 


u t + Ux = o X £ (0, 1) 
u(z,0) = f(x) 
u{0,t) = g(t) 

The corresponding semi-discrete system reads 


v t + PDv - {I - P)g t 

v(0) = / 


0 

\ 


9 

\ 

1 



9 1 


• 

. 

g = 




1 J 


9u 

/ 


(71) 


(72) 
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where gj, j = are to be determined later. Obviously, Vo (t) = g(t) if / 0 = (/(O) 

(cf. proposition 2.4). The boundary condition is thus fulfilled at all time. According to 
proposition 2.4 one has v = Pv + (I - P)g (/ - P)(v - g) = 0. Hence, the energy 
method gives 

j t \Hl = -2(v -(I- P)g, Dv) h + 2(v,(/ - P)g,) h 
Subtracting 2 (g, v t )h from both sides we get 


2(v - g, v t ) h = —2(v - g, Dv) h - 2 (g, v t + PDv - {I - P)g t )h. + 2(u - g,(I - P)g t )h 


Using (72) and (I — P)(v — g) = 0 shows that 


2(v - g , v t )k = —2(v - g , Dv) h 


i. e., 

2(v -g,(v- g)t)h = -2(u - g , D(v - g)) h - 2(v - g,g t + Dg) h 
If q solves the auxiliary problem 

9t + Dg = 0 

9 ( 0 ) = / 

then 

^\\v -g\\l = -2(v - g, D(v - g)) h = (v 0 -gf - K ~ 9v? < 0 
since = 5 . Thus, 

\\v(t)-m\\H<\\v(o)-m\\H = o 


Consequently, 

v(t) = g(t ), 

If g satisfies (73) we get the energy estimate 


t > 0 


j t M\l = -2 (g,Dg)h = g 2 -gl 

Hence, 

\\m\\U fgl(T)dT= ||/||fc+ fg\r)dr 
Jo Jo 

Finally, v = g implies 

IM0IG+ /'i>;w*- = ||/||i+ fg 2 (T)dr 

Jo Jo 


(73) 


which is identical to the continuous estimate. 

It remains to be verified under what circumstances g solves the auxiliary problem (73). 
Let g — ( go g\ ... g v ) T be the solution to (73). Hence 


D’g(0) = D 3 f j = 0,1,-.. 
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and 


+ (-l) j+1 D j g = 0, i>0, j = 0, 1,... 


i.e., for t = 0 we get the compatibility conditions 


+ = 0, j = 0,1,... 

Thus, if we require that the initial-boundary data / and g satisfy 


•(0) + (-i) j+1 (£>V)o = 0, i = o, l, . . . 


it follows that 

dV ^ dV ^ 3 ’ ’ ‘ ‘ ' 

By virtue of being the solution to (73) go(t) is analytic in t. Hence, taking the boundary 
data g(t) to be analytic these equalities imply that g(t) = ^ > 0, i.e., g = g, which 

proves that g indeed solves (73). 

In what follows we shall analyze the general case. Consider the ODE-system 

(J-'v), = PR(t,v) + (I-P)(J-'g), 

v(0 ) = / ( ' 4) 

with J _1 being the inverse Jacobian, and where 

R(t, v) = G(t, v) + r l F{t ), G(t, 0) = 0 

This form arises naturally when discretizing a non-linear PDE in space; g represents the 
boundary data, and F is the forcing function; G(t,v) is the discretization of the differen- 
tial operator. It should be pointed out that most operators G occurring in practice are 
autonomous , i. e., G = G(v). We use the tilde notation to emphasize that g is only par- 
tially determined, that is, some components are determined by the boundary data of the 
underlying PDE; the remaining components are unknown. It is no restriction to assume 
that g = ( go ... g u ) T with i = 0, . . . , s being the known components. Otherwise, 
g could be brought to this form by permuting the dependent variables appropriately. As 
usual we require that P and J~ l commute, which is true if the grid is locally isochoric at 
the boundary. Next, we define the auxiliary problem 

(J~ l w) t = R(t,w) / 75 ) 

u>(0) = / { } 


Any solution to (75) will satisfy 


(j 1 w) = Rj(t,w), j = 0,1,... 
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where R } is defined recursively by 

f)R BR 

Rj{t,w) = ( t,w)R(t,w ), j = 1,2,... 

Ro(t, w) = w 

Consequently, at t = 0 we have 

= j = 0,1,... 

Assumption 6.1 The boundary data gi(t), i = 0, . . . , s, the initial data f , and the forcing 
function F satisfy the compatibility conditions 

(j~ (0) = (^(0, /)),-, j = 0,...,5, j = 0, 1, . . . 

If R is sufficiently well-behaved, in particular if G is linear and autonomous, then w(t) will 
be analytic for 0 < t < T . Thus, if we require that the boundary data gift), i — 0, . . . , s, 
be analytic it follows that 

gift) = Wi(t), i = 0, ...s 

Furthermore, the unknown components i > s are of course taken to be 

g t (t) = Wi(t ), i > s 


Hence, g = w solves (75). 

Remark: It suffices to consider g x piecewise analytic, since the process can be repeated 
at t = where G is the critical time when analyticity is lost. 


Proposition 6.1 Let v be the solution to (74) and suppose that assumption 6 A holds. If 
the boundary data g{ , i = 0, . . . , s are piecewise analytic, then 

{v - g,{v - g)t) h = (t> - g,G{t,v) - G(t,g)) h 


Proof: 

Using (74) and Pv = v — (I — P)g , which is true since P and commute, it follows 
readily that 


(v, v t )h = (v-g, R(i, v)) h + (< 7 , PR(t , v)) h + (u, (/ - P){J 1 g)t)h 


Hence, 


(v-g,v t )h = {v-g, R(t,v)) h + (g,-{J 1 u) ( + PR(t, v)) h + (w, (I - P){J 1 g) t )h 
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or, using (74), 


(v-g,v t ) h = (v - g, R(t,v)) h + (v - g,{I - P)(J 1 g) t )h 

But 

(v - h (/ - P){J~ x 9 )t)h = ((/ - P)(v - g), (J-"g)t)k = 0 

and so 

{v - g, v t )h = {v-g, R(t, v)) h 

which in turn is equivalent to 

(v - g, ( v-g)t)h = (v - g,R(t,v) - R{t,g)) h ~ (w - g)t - R(t,g))h 
The assumptions on g imply that (J~ x g) t = R(t,g), which proves the proposition. □ 

6.1 One-Dimensional Parabolic Systems 

We consider the parabolic system (27) with the lower order terms omitted, i. e., 

lit — (^d^x)a: T P 
u{x,0) = f 

Lou(0,t) + Liu x (0,t) = g(t) 

The omission of lower order terms is done for convenience only. The boundary data g(t) 
is assumed to be piecewise analytic in t. The corresponding semi-discrete system reads 

v t = P ( DADv + F) + (I- P)g t , . 

v(0) = f 

where g satisfies L T g = g , with L given by (18). According to proposition 2.4 we have 
(/ — P)(v — g) or, equivalently, L T v = L T g. Thus 

L 0 vq + L](Dv) 0 = g 

which shows that the analytic boundary conditions are satisfied to some order of accuracy. 
Proposition 6.2 Suppose that assumption 6.1 holds. Then (77) is strictly stable. 

Proof: 

We know that g solves the auxiliary problem 

9t = £{*,</) + F (t) 

m = / 


m = ( 9 o*’ ) < 76) 
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where G(t,g ) = DADg. Proposition 6.1 then yields (using J = I) 

(v-g, {v-g) t )h = (v - g,DAD(v - g)) h < ~{v 0 - g 0 ) T A(D{v - g)) Q - 28\\D{v - g)\\\ 
Since L T (v — g) — 0 it follows that 

I 0 >0 - £o) = o 

Furthermore, decompose vj - g 3 = (u' — g'-) + (v" - g"), where u' , g' } € ker(L{), u", g" € 
ker^l) 1 . According to assumption 3.1 we then obtain 

(v - g, (u - 5)t)fc < -(«b - go? A(D{v" - g")) 0 - 2£||£>(t7 - <7)|| 2 

Arguing exactly as in the proof of proposition 3.2 gives 

(D(v"-g ")) 0 = -LT 1 Lo(v 0 -go) 


i- e., 

(v - g,(v - g) t ) h < ilvo - g 0 \ 2 - 28\\D{v - g)\\l 
By means of the Sobolev inequality 2.6 we thus arrive at 

(v -g,(v- g)t)h < (<y + 0{h))\\v - g\\l 

Hence 

IMO - 5(011* < e <“+°<'»>'||t,( 0 ) - j( 0 )|U = o 

which is equivalent to v(t) = t > 0. To get the final estimate we consider the 

auxiliary problem. One obtains 

j t M\l < -2g r 0 A(Dg")o ~ 4<||Dj|i; + + \\F\\l 

Logo + Li(Dg")o = g 
(Dg")o = -i^ 1 L 0 go + L~ l g 


Now 
and so 
Thus 


-2goA(Dg') 0 = 2 gl A~L X 1 L 0 g 0 - 2gl AL j l g < 7|p 0 | 2 + |tf| 2 
where the algebraic inequality 2xy < ex 2 + e~ 1 y 2 was used. This leads to 

J t m\l + M* < (7 + l)l«o| ! - ^IIDslli + llsllJ + Iff! 2 + IICIII 

The coefficients of this estimate are exactly the same as those of the corresponding analytic 
inequality. Using g = v and eliminating the boundary terms of the right member by means 
of the Sobolev inequality yields 

^IMI 2 + M 2 < (« + £ > (M)IM| 2 + \g\ 2 + HF’II 2 
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where a is the same constant as in the analytic estimate. Finally, integration with respect 
to time results in 

Ml + j ‘ Mr)| 2 dr < (||/||J + £ (| S (r)p + ||F(r)||J) dr) 

which is the desired estimate. □ 

Remark: The boundary conditions are used twice — first in conjunction with propo- 
sition 6.1 to show that g = v, and second with the auxiliary problem to get the actual 
estimate. 

6.2 Two-Dimensional Symmetric Hyperbolic Systems 

We consider (the lower order terms are omitted for convenience) 

(j-'u), = \ £ + J~'F (78) 

where 

Bi = dti • A = pi- At + |~ A 2 

UXi OX 2 

The boundary conditions are given by 

<pi(x,t) = S(x)<pn(x,t) + g(x,t) 

where g>\ designates the characteristic variables corresponding to the locally ingoing char- 
acteristics; S(x) is assumed to be sufficiently small. At the corner x(0,0) we require 
that 

y>/(x(0,0),f) = S(x(0,0))v?//(x(0,0),f) + y(x(0,0),f), ^>(x(0, 0), t) = Qju(x(0,0),t) 

be satisfied for i = 1,2, where C2f(rci^-Ai(x(0, 0),t) + ra 2 ^.A 2 (x( 0 , 0),f)Q,- are diagonal; 
i = 1,2 are the two normals associated with the corner x(0, 0). 

At each discrete boundary point x±j = x{ih\^jh£) G T the boundary conditions are 
formulated as 

T(x,j )u,j = <7(xf,,t), X{j G T 

where L(x,y) is given by (32). We note that there are two operators Li(x 0 0 ) corresponding 
to the two different normals n b) at xoo- Since L(x tJ ) has full rank it follows that there 
exists a g tJ such that 

L(xij)gij = g(xij,t), Xij g T 

Hence, the boundary conditions are 

T( X,j ) ( Vjj 9ij) ^ r 
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or, in global form, 


L T (v -g)= 0 

where g is partially determined by the boundary data. Eq. (78) is discretized as 

(j-'v), = p (5 £ (D,J-'B;v + J-'B,D,v) +J-'F)+(1-P) (. r'g \ (7g) 

v(0) = 1 

where P is the orthogonal projection corresponding to the global operator L r . 
Proposition 6.3 Suppose that assumption 6.1 holds. Then (79) is strictly stable. 

Proof: 

By assumption 6.1 g solves the auxiliary problem 
where 

G(t,g) = \t (DiJ-'B# + J-'BiDtg ) 

1 1 = 1 

Hence, according to proposition 6.1 we have 

(v - g,{v- g)t)h = \ J2( v - 9, (DiJ~ l Bi + J-'BiDi) (v - g)) h 

Summation by parts yields 

(v -g,(v- g)t)h = l(u - 5 , (njAi + n 2 A 2 ) {v - g))r 

But L T (v — g) = 0, whence v — g satisfies the homogeneous boundary conditions. Conse- 
quently, (cf. the proofs of propositions 3.1, 4.1, 5.1) 

(v - g , (nji4i + n 2 A 2 ) (v - g)) r < ~^\v - < 0 

Thus, g(t ) = v(t), t > 0. 

In the second part of the proof we apply the energy method to the auxiliary problem. 
Straightforward computations show that 

4 (9,g)h = {g,{nxAi + n 2 A 2 ) g) r + 2{g, F) h 

Take an arbitrary point x 0 j on the boundary portion where = 0. We must analyze the 
quadratic form 

9oj + n^A^.go, 
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We know that L T g — g, where g is the vector representing the analytic boundary data. 
Define <p tJ — Q T (x,j)gij. Hence, 

= S(xij)(¥ij)u + g tJ (80) 

and the quadratic form is transformed into 

Soj ( n i 1} ^i + A 2 ) fa = Vojh-onpoj 

Using (80) gives (omitting the spatial subscripts for simplicity) 

ip T Ag? = <fpjj (A u + S t AjS^ (fn + 2<pJjS T Aig + g T Ajg 
It is assumed that A// < —7 at x Gj . For sufficiently small S we thus get 

< -jlw/l 2 + (l + |A,|) |j | 2 

Now, |(^/| < |5||y>//| + Hence 

<? T A<p + ^|^| 2 < -jlvn] 2 + ^l^/| 2 + (1 + l A /|) \g\ 2 < (3 + | A/|) \g\ 2 

for S small enough. It should be underscored that this is exactly the same estimate one 
gets in the continuous case. At each boundary point Xo j we have thus established that 

9oj + n ( 2 ] A 2 ) o .g 0j + ^p|<7oj| 2 < ( 3 + l( A oj)/|) Itfojl 2 

with a similar expression at points x,q corresponding to £2 = 0. Letting inf( 7 ,j) = 7 > 0 
we thus obtain 

j t {g,g)h + ^{g,g ) r < (3 + \\\oo){g,g)r + 2{g,F) h 
Finally, identifying v = g, integration gives the energy estimate 

(»(<),«(<))» + jfV(r), v(t)) t , dr < Ke‘ ((/, f) k + J‘ «j(t), j(t)) p + (F(r), F(r))») dr) 

which proves the theorem. D 

Remark: Because of the terms <pfjS T Aig + g T Aig the constant K of the energy estimate 
will in general satisfy K > 1 , even if no estimate of the boundary terms (v,v)r is wanted. 
For g = 0, i. e., homogeneous boundary conditions, the critical terms disappear, and one 
may take K = 1 in case no boundary estimate is needed (cf. the remark following the 
proof of proposition 3.1). 
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7 Summary and Conclusions 


We have demonstrated that for a given finite-dimensional scalar product (•, -)h any linear 
discretized boundary condition can be written as an orthogonal projection operator P that 
satisfies (u, Pv)h = ( Pu,v)h . It should be noted that the projection is well-defined if the 
corresponding analytic problem is well-posed. For general boundary conditions one may 
also have to require that the discretization parameter h be small enough (consistency). 
The projections P , the summation-by-parts property, and proposition 2.4 constitute the 
main tools needed to obtain an energy estimate for the semi-discrete case. For a large 
class of problems it has been established that existence of an energy estimate for the 
continuous problem implies ditto for the semi-discrete system. 

In one space dimension we are no longer required to consider restricted full norms 



which were used in [4] to prove stability for symmetric hyperbolic systems subject to 
homogeneous boundary conditions. The main result is the stability proof for mixed 
hyperbolic-parabolic systems subject to general linear boundary conditions. Assuming 
certain compatibility conditions the result holds for inhomogeneous boundary data. Re- 
formulating the analytic problem it is possible to obtain strict stability, i. e., we have 
a time stable semi-discrete approximation that is bounded by the the same exponential 
growth rate (modulo O(h)) as the analytic problem. For the parabolic part the excess 
growth rate is induced by the discrete Sobolev inequality. Furthermore, for the hyperbolic 
part we have used assumption 2.1. In particular, strict stability is obtained for diagonal 
norms and variable coefficient problems, and for general norms and constant coefficient 
problems. The stability results hold for finite difference approximations of arbitrary order. 

In two space dimensions we are forced to consider diagonal norms in order to get 
summation by parts in both dimensions. Stability of high-order schemes is obtained for 
general mixed hyperbolic-parabolic initial-boundary value problems. Again, inhomoge- 
neous boundary conditions are allowed, provided certain compatibility conditions prevail. 
Using a different norm we obtain strict stability for symmetric hyperbolic systems on 
non-smooth curvilinear domains, where we allow for general inhomogeneous boundary 
conditions. As for strict stability of parabolic systems, we are limited to homogeneous 
Dirichlet data. Mixed derivatives and/or variable geometry may account for “tangential 
differences” that in general cannot be eliminated without ruining strict stability. An ex- 
ception is the standard second order method. The requirement that the Dirichlet data 
be homogeneous is already necessary for the continuous problem. All results obtained for 
two dimensions generalize to higher dimensions. 

The methods presented in this report are similar to finite element methods in that 
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stability for the semi-discrete system follows more or less directly from the corresponding 
continuous one. There is, however, one major difference: The FEM technique often results 
in implicit space discretization, whereas the discretized space operators reported in this 
paper always are explicit. 

There are other ways of imposing boundary conditions so as to ensure time stability 
(strict stability) when using difference operators satisfying a summation-by-parts prop- 
erty. An elegant technique is proposed in [1]. A so called Simultaneous Approximation 
Term, SAT for short, is added to the semi-discrete scheme. The SAT will act as a penalty 
function to enforce an approximation of the discrete boundary conditions. In [1] this 
approach is used to prove time stability for high-order finite difference approximations 
of one-dimensional constant coefficient hyperbolic systems. Also, it is not necessary to 
consider identical difference stencils in the interior. A new and interesting class of such 
difference operators can be found in [2]. 
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8 Appendix 


Let P be defined by proposition 2.5 with L given by (42) (x = 0.5) where h\ = h ,2 = \ 
for convenience. We shall show that Lj 0 P = 0, which will follow if we can prove that 
L(L t T ,~ 1 L)~ r L t = Lj 0 . Straightforward computations show that 
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Obviously, 1 / is the Schur complement of D 3 2 - Let L be defined by (40). Then 
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where (using T X2 = T 2 3 and T u = T 33 ) 


( T 22 /A T?J2 \ ( T 22 /4 Tj 2 / 2 \ 

l T u / 2 T„ j ^ T„/2 T i3 J 


Furthermore, 

L[ 0 S" 1 L= ( Ac/ao 0 4 /(t 0 2 2£>f 2 ) 

Using D 22 = (k/ctq + dl 0 /cr%)/ 2 , 
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d 22 t? 2 + 2£»f 2 r 13 = 2 (Df 2 T 13 + D 22 Tl + ^ 2 r„) - + 2Df a r„) 

Observing that 7^ = T 23 , Dj 2 T u = D 23 T 33 it follows that the first parenthetical expres- 
sion vanishes. Thus, 
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Also, 
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Substituting 1 — err 2 = dl 0 /(Ka 0 ) and the expression for /x yields 
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Consequently, 
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The j th block column of L T reads 
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where the sum is the j : th block element of Cj, j = 0, . . . ,r. Hence, 


r-j — ( <^00 —doi ■ ■ ■ —dor doo doi ■ ■ ■ do r ) c j — ^?o ^ d ok e k 
Zd 0 o k=0 

Lf 0 S- 1 L(I r S“ 1 L)- 1 L T = Lj, 

In a similar fashion one shows that L 20 also satisfies the above equation. 

The simplest example is obtained by discretizing the Neumann conditions using the 
standard divided difference D+ in both coordinate directions, i. e., 

^ii — v oi = 0 

t>n-uio = 0 (81) 

(uio — voo)/2 + (voi _ vqo)/2 = 0 

which leads to 

16 -4 0 ... 0 -4 -8 0 ... 0 

-2 14 0 ... 0 -4 -8 0 ... 0 

0 0 0 ... 0 0 0 0 ... 0 

0 0 0 ... 0 0 0 0 ... 0 

-2 -4 0 ... 0 14 -8 0 ... 0 

-2 -4 0 ... 0 -4 10 0 ... 0 

0 0 0 ... 0 0 0 0 ... 0 

0 0 0 ... 0 0 0 0 ... 0 

Evidently, any vector v = Pu, P = I — E -1 L(L r E _1 L) _1 Z/ T , satisfies (81). Furthermore, 
by (82), V 10 — vqo = ^01 — Voo = 0, which also follows directly from (81). 
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